213-Discussion

Discussion I

For today’s discussion, we will review the distinction between parametric and nonparametric regression models, as well as some key properties of simple linear regression. If you wish to review materials from POL 211 & 212, you can access them through the following links: POL 211 Discussion, POL 212 Discussion.

1. Parametric versus nonparametric Regression Models
2. Simple Linear Regression
3. Properties of Least Squares Estimator

1. Parametric versus nonparametric Regression Models

When it comes to regression analysis, choosing the right approach is crucial for accurate predictions and meaningful insights. Two common methods used are parametric, like linear regression, and semi/non-parametric, like smoothing spline regression or Kernal regression. Each has its own advantages and disadvantages, and the choice between them largely depends on the nature of the data and the underlying relationships.

a. What is parametric and semi/nonparametric regression?

Parametric Regression Linear regression is a well-known parametric method that assumes a linear functional form for the relationship between the predictors (\(X\)) and the target variable (\(Y\)). This approach has several benefits, such as ease of estimation with a small number of coefficients. In linear regression, these coefficients have straightforward interpretations, and statistical significance tests are readily applicable. However, parametric methods come with a significant limitation — they rely on the assumption that the specified functional form is a close approximation to the true relationship. If this assumption is far from reality, linear regression can perform poorly and yield unreliable results.

Nonparametric Regression On the other hand, non-parametric methods like K-Nearest Neighbors (KNN) regression do not make explicit assumptions about the functional form of the relationship between \(X\) and \(Y\). Instead, they provide a more flexible approach for regression. KNN regression identifies the K training observations closest to a prediction point and estimates the target variable by averaging their values. While this approach is more versatile and can handle complex relationships, it can suffer from high variance when K is small, leading to overfitting. Conversely, when K is large, KNN regression can underfit the data.

Example

Assume that we have a outcome variable \(Y\) and two explanatory variables, \(x_1\) and \(x_2\). In general, the regression model that describes the relationship can be written as:

\[Y = f_1(x_1) + f_2(x_2) + \epsilon\]

• Some parametric regression models:

  • Multiple linear regression model: \(Y = \beta_0 + \beta_1x_1 + \beta_2x_2 + \epsilon\)
  • Polynomial regression model of second order: \(Y = \beta_0 + \beta_{10}x_1 + \beta_{11}{x_1}^2 + \beta_{20}x_2 + \beta_{21}{x_2}^2 + \epsilon\)
  • Nonlinear regression model: \(Y = \beta_0 + \beta_1x_1 + \beta_2e^{\beta_3x_2} + \epsilon\)
  • Poisson regression with \(Y\) is count: \(log(\mu) = \beta_0 + \beta_1x_1 + \beta_2x_2 + \epsilon\)

• If we do not know \(f_1\) and \(f_2\) functions, we need to use a Nonparametric regression model.

b. Nonparametric regression estimation methods: K-Nearest Neighbors (KNN) regression

K-Nearest Neighbors (KNN) regression is one of the simplest and best-known nonparametric methods. Given a value for \(K\) and a prediction point \(x_0\), KNN regression first identifies the \(K\) training observations that are closest to \(x_0\), represented by \(N_0\). It then estimates \(f(x_0)\) using the average of all the training responses in \(N_0\). In other words,

\[\hat{f}(x_0)=\frac{1}{K}\sum_{x_i \in N_0}y_i\]

When using KNN as a regressor with a continuous dependent variable, data points are scattered across the coordinate plane. When a new data point is introduced, the number of neighbors (\(K\)) is determined using any of the distance metrics.Usually, the Euclidean distance is used as the distance metric. Once the neighbors are identified, the predicted value of the new data point is calculated as the average of all the neighbors’ values combined.

The below figure illustrates two KNN fits on a data set with \(p = 2\) predictors. The fit with \(K = 1\) is shown in the left-hand panel, while the right-hand panel corresponds to \(K = 9\).We see that when \(K = 1\), the KNN fit perfectly interpolates the training observations, and consequently takes the form of a step function. When \(K = 9\), the KNN fit still is a step function, but averaging over nine observations results in much smaller regions of constant prediction, and consequently a smoother fit. In general, the optimal value for \(K\) will depend on the bias-variance tradeoff, which Chris introduced in POL 212. A small value for \(K\) provides the most flexible fit, which will have low bias but high variance. This variance is due to the fact that the prediction in a given region is entirely dependent on just one observation. In contrast, larger values of \(K\) provide a smoother and less variable fit; the prediction in a region is an average of several points, and so changing one observation has a smaller effect. However, the smoothing may cause bias by masking some of the structure in \(f(X)\).

The bias-variance tradeoff is a fundamental concept in machine learning that deals with the balance between the bias (error from overly simplistic assumptions) and variance (sensitivity to fluctuations in the training data) of a model. In the context of KNN (K-Nearest Neighbors), adjusting the value of K can help control this tradeoff.

• Low \(K\) (e.g., \(K=1\)): Low bias, high variance. The model closely follows the training data, leading to high variance and sensitivity to noise. It can capture complex patterns but may overfit the training data.

• High \(K\) (e.g., large \(K\)): High bias, low variance. The model averages over more data points, resulting in smoother predictions with lower variance but potentially higher bias. It may underfit the data by oversimplifying the underlying patterns.

Choosing an appropriate value for \(K\) involves finding a balance between bias and variance to achieve optimal model performance. This tradeoff influences the model’s ability to generalize to unseen data.

There are no pre-defined statistical methods to find the most favorable value of \(K\), but here is a typical way of how to choose an optimal \(K\):

1. Initialize a random \(K\) value and start computing.
2. Choosing a small value of \(K\) leads to unstable decision boundaries.
3. The substantial \(K\) value is better for classification as it leads to smoothening the decision boundaries.
4. Derive a plot between error rate and \(K\) denoting values in a defined range. Then choose the \(K\) value as having a minimum error rate.

c. When should we use parametric or nonparametric regression?

The key question is when to choose a parametric approach like linear regression over a non-parametric one such as KNN regression. The answer is straightforward: a parametric approach performs better when the chosen functional form is a close match to the true relationship, particularly in the presence of a linear relationship. If the specified functional form is far from the truth, and prediction accuracy is our goal, then the parametric method will perform poorly. For instance, if we assume a linear relationship between \(X\) and \(Y\) but the true relationship is far from linear, then the resulting model will provide a poor fit to the data, and any conclusions drawn from it will be suspect.

In contrast, non-parametric methods do not explicitly assume a parametric form for \(f(X)\), and thereby provide an alternative and more flexible approach for performing regression.

To illustrate this point, let’s consider a few scenarios and use R to simulate them:

1. Linear Relationship: When the true relationship between \(X\) and \(Y\) is linear, linear regression outperforms nonparametric regression. Linear regression provides an almost perfect fit in this situation, as it closely matches the underlying relationship.

   set.seed(1234) # for replication
   N <- 1000 # set up the sample size
   x1 <- rnorm(N, mean = 10, sd = 3)
   e <- rnorm(N, mean = 0, sd = 3) # set the errors to be normally distributed 
   y <- 7 + 3*x1 + e # set the true y
   df <- data.frame(id = 1:N, y = y, x1 = x1, e = e)

   In the above chunk, I specify the true \(y\) to be a function of \(x_1\) and \(e\). The true values of the regression coefficients is 3 for \(x_1\). The intercept is 7. Let’s first plot the true relationship between \(x_1\) and \(y\)!

   library(tidyverse)
   ggplot(data = df, aes(y = y, x = x1)) +
     geom_point(alpha = 0.6) +
     geom_smooth() +
     theme_bw()

   

   Let’s begin by splitting the generated data into a training set, comprising 80% of the data, and a test set, comprising the remaining 20%.

   set.seed(1234) # set the seed to make the partition reproducible
   train <- df %>% sample_frac(0.8) # should have 800 obs
   test <- anti_join(df, train, by = 'id') # should have 200 obs

   Now, let’s first apply the lm() function to our training set. We’ll then use the coefficients estimated from this model to predict the y-values in the test set, followed by computing the mean squared error (MSE).

   fit <- lm(y ~ x1, data = train)

   	OLS
   Variables	Estimates	S.E.	C.I. (95%)
   (Intercept)	6.50 ***	0.35	5.81 – 7.19
   x1	3.06 ***	0.03	2.99 – 3.13
   Observations	800
   R2 / R2 adjusted	0.910 / 0.910
   p<0.05   ** p<0.01   *** p<0.001

   test$y_hat <- predict(fit, newdata = data.frame(x1 = test$x1))
   
   # compute mean squared error (MSE)
   MSE_ols <- sum((test$y - test$y_hat)^2)/nrow(test)
   
   MSE_ols

   ## [1] 8.770813

   The mean squared error (MSE) obtained by using OLS regression to estimate the relationship between \(y\) and \(x_1\) in the training set and then using the estimated coefficients to predict \(y\) in the test set is 8.770813. Now, let’s proceed to use KNN regression with various values of \(K\) to observe how the MSE varies.

   We will use the knn.reg() function from the FNN package for KNN regression. This function requires four arguments to be specified:

   • train: the predictors of the training data
   • test: the predictor values, \(x\), at which we would like to make predictions
   • y: the response for the training data
   • k: the number of neighbors to consider

   # install.packages("FNN")
   library(FNN)
   train_no_y <- train %>% select(x1) 
   test_no_y <- test %>% select(x1) 
   train_y <- train %>% select(y) 
   
   yhat_k1 <- knn.reg(train = train_no_y, test = test_no_y, y = train_y, k = 1)
   yhat_k1

   ## Prediction:
   ##   [1] 29.381821 28.304426 65.609505 23.837209 34.330518 33.436274 21.824014
   ##   [8] 30.847990 32.289959 42.276634 48.391432 53.703199 31.595459 33.423440
   ##  [15] 32.935792 41.304789 26.486024 39.445252 52.119807 38.608137 44.000824
   ##  [22] 38.608137 37.159409 38.540173 33.079801 39.907787 44.810978 37.870818
   ##  [29] 34.348167 42.174368 40.401263 34.996959 27.478094 48.391432 33.978417
   ##  [36] 26.367731 56.238177 31.240796 37.950729 55.185947 34.850824 43.470660
   ##  [43] 44.483329 21.700094 49.108239 44.558484 41.288458 40.623165 31.965684
   ##  [50] 28.174148 35.309716 44.393660 40.512042 48.391432 26.800229 31.877765
   ##  [57] 45.551430 28.354259 39.391663 43.100323 33.168206 25.652004 37.857259
   ##  [64] 42.280182 31.905242 30.070676 45.434786 32.597936 55.713238 30.245094
   ##  [71] 35.747199 24.318503 35.309716 42.142535 44.687740 34.463743 41.607736
   ##  [78] 11.887226 36.074342 55.185947 33.106591 64.709783 20.566998 54.842595
   ##  [85] 45.434786 51.245243 39.073328 28.874519 41.177921 35.357025 30.071689
   ##  [92] 36.168858 35.095983 28.392045 31.241205 32.798233 29.709591 37.191013
   ##  [99] 38.608137 26.841422 44.442204 24.630053 26.367731 19.105273 33.290191
   ## [106] 40.401263 56.238177 21.708475 33.771862 38.654251 42.868435 49.161958
   ## [113] 33.168206 26.492572 27.573812 41.597145 30.479421 45.411737 40.623165
   ## [120] 37.208831 35.170920  8.539242 48.085100 44.483329 20.591604 21.824014
   ## [127] 45.068575 40.911258 35.964535 44.731334 38.432659 44.912912 38.520474
   ## [134] 33.390623 24.145923 22.517012 37.950729 41.465057 36.782189 31.241205
   ## [141] 33.423440 29.709591 28.827443 21.824014 39.739236 49.108239 41.337496
   ## [148] 51.204406 31.801673 31.406837 30.414360 32.690670 39.641522 28.351906
   ## [155] 44.206919 32.373809 48.402430 31.595459 67.724169 27.404886 36.782189
   ## [162] 42.020108 54.160505 36.074342 48.649611 34.330518 44.253326 42.549677
   ## [169] 36.377909 46.279074 27.157538 25.692918 30.150333 29.936407 50.422834
   ## [176] 27.404886 37.870818 33.978417 29.820616 32.646670 28.464324 26.137602
   ## [183] 29.774837 38.787417 29.728468 39.795625 54.842595 48.649611 38.169491
   ## [190] 38.146542 33.233128 36.284027 41.820483 35.838288 38.441584 32.682236
   ## [197] 28.521211 41.054304 52.592726 44.479609

   As observed from the above results, when setting \(K = 1\) in knn.reg(), it provides the predicted \(y\) values for the test set. Now, let’s calculate the MSE for this prediction.

   MSE_k1 <- sum((test$y - yhat_k1$pred)^2)/nrow(test)
   MSE_k1

   ## [1] 20.25002

   As you can see, the MSE when employing KNN regression with \(K = 1\) is significantly higher than the MSE observed with OLS regression. Let’s create a for loop and repeat this process for various values of \(K\) to determine whether the MSE improves or worsens compared to the MSE when using OLS.

   MSE_KNN <- c()
   
   for (i in seq(1, 100, by = 2)){
     yhat_k <- knn.reg(train = train_no_y, test = test_no_y, y = train_y, k = i)
     MSE_k <- sum((test$y - yhat_k$pred)^2)/nrow(test)
     MSE_KNN <- c(MSE_KNN, MSE_k)
   }
   
   MSE_KNN

   ##  [1] 20.250016 11.554227 10.170821 10.020766 10.002982  9.655185  9.468217
   ##  [8]  9.151752  9.058586  8.825635  8.704344  8.733093  8.672769  8.708741
   ## [15]  8.678995  8.627246  8.531880  8.596169  8.653730  8.691982  8.783636
   ## [22]  8.831042  8.858023  8.874769  8.833497  8.853275  8.927840  8.982527
   ## [29]  8.938121  8.948997  9.004996  9.012363  8.970862  8.969159  9.031482
   ## [36]  9.058308  9.110250  9.158928  9.255961  9.305864  9.352528  9.415657
   ## [43]  9.468395  9.539323  9.576483  9.650497  9.693073  9.724640  9.752314
   ## [50]  9.824734

   ggplot(data = data.frame(K = seq(1, 100, by = 2), MSE = MSE_KNN), 
      aes(y = MSE, x = K)) +
     geom_line(alpha = 0.4, color = "darkgreen") +
     geom_point(alpha = 0.3, color = "darkgreen") +
     geom_hline(yintercept = MSE_ols, color = "darkgray", 
            alpha = 0.8, linetype = "dashed") +
     annotate("text", x = 90, y = 18, label = "KNN", color = "darkgreen", alpha = 0.5) +
     annotate("text", x = 90, y = 17, label = "OLS", color = "darkgray") +
     theme_bw()

   

   After simulation, it’s evident that when the true relationship between the variables of interest is linear, KNN regression doesn’t offer superior predictions compared to OLS regression.

2. Slight Non-Linearity: In cases of slight non-linearity, where the true relationship deviates slightly from linearity, KNN regression can perform nearly as well as linear regression. It still provides reasonable results without a substantial reduction in prediction accuracy.

3. Strong Non-Linearity: However, in situations with a strong non-linear relationship, KNN regression outperforms linear regression. This is because KNN can adapt to complex relationships, providing more accurate predictions.

   set.seed(1234) # for replication
   N <- 1000 # set up the sample size
   x1 <- rnorm(N, mean = 1, sd = 10)
   e <- rnorm(N, mean = 3, sd = 10) # set the errors to be normally distributed 
   y <- 3 + 4*(x1+3)^2*(x1-1)^3 + e # set the true y
   df <- data.frame(id = 1:N, y = y, x1 = x1, e = e)

   For example, let’s now generate a relationship between \(x_1\) and \(y\) that is strongly non-linear. I specify the true \(y\) to be a function of \(x_1\) and \(e\). Overall, the formula, \(y = 3+4(x_1+3)^2(x_1-1)^3+e\), describes a relationship where \(y\) is influenced by \(x_1\) in a nonlinear way, with two distinct “bumps” or curves, along with some random variability represented by \(e\).

   library(tidyverse)
   ggplot(data = df, aes(y = y, x = x1)) +
     geom_point(alpha = 0.6) +
     geom_smooth() +
     theme_bw()

   

   Following our previous approach, let’s now compare the MSE obtained from OLS regression, used to estimate the relationship between \(y\) and \(x_1\), with the MSE derived from KNN regression, employing different values of \(K\).

   • Step 1: Split the Sample in to Training and Test Set

     set.seed(1234) # set the seed to make the partition reproducible
     train <- df %>% sample_frac(0.8) # should have 800 obs
     test <- anti_join(df, train, by = 'id') # should have 200 obs
     fit2 <- lm(y ~ x1, data = train)

   • Step 2: Run OLS Regression and Compute the MSE

     fit2 <- lm(y ~ x1, data = train)

     	OLS
     Variables	Estimates	S.E.	C.I. (95%)
     (Intercept)	95970.43	389692.62	-668973.27 – 860914.13
     x1	734839.02 ***	38319.47	659620.15 – 810057.88
     Observations	800
     R2 / R2 adjusted	0.315 / 0.315
     p<0.05   ** p<0.01   *** p<0.001

     test$y_hat <- predict(fit, newdata = data.frame(x1 = test$x1))
     MSE_ols <- sum((test$y - test$y_hat)^2)/nrow(test)
     
     MSE_ols

     ## [1] 4.742095e+13

   • Step 3: Run KNN and Compute the MSE

     train_no_y <- train %>% select(x1) 
     test_no_y <- test %>% select(x1) 
     train_y <- train %>% select(y) 

     MSE_KNN <- c()
     
     for (i in seq(1, 100, by = 2)){
       yhat_k <- knn.reg(train = train_no_y, test = test_no_y, y = train_y, k = i)
       MSE_k <- sum((test$y - yhat_k$pred)^2)/nrow(test)
       MSE_KNN <- c(MSE_KNN, MSE_k)
     }
     
     MSE_KNN

     ##  [1] 1.006088e+11 1.791096e+11 1.840990e+11 5.863832e+11 9.814868e+11
     ##  [6] 4.548076e+11 1.081960e+12 1.129925e+12 1.487238e+12 2.374122e+12
     ## [11] 2.763030e+12 3.688595e+12 4.638701e+12 5.402272e+12 6.302316e+12
     ## [16] 6.098534e+12 6.924194e+12 7.720989e+12 8.508775e+12 9.279200e+12
     ## [21] 9.522820e+12 1.021639e+13 1.021119e+13 1.086252e+13 1.149660e+13
     ## [26] 1.209597e+13 1.238768e+13 1.296902e+13 1.354249e+13 1.352895e+13
     ## [31] 1.408380e+13 1.458944e+13 1.475359e+13 1.524066e+13 1.571534e+13
     ## [36] 1.618749e+13 1.664102e+13 1.702627e+13 1.745704e+13 1.775741e+13
     ## [41] 1.815320e+13 1.818130e+13 1.853354e+13 1.892895e+13 1.930142e+13
     ## [46] 1.966064e+13 2.001759e+13 2.037685e+13 2.072546e+13 2.106549e+13

   • Step 4: Compare the MSE obtained from OLS and with that from KNN

     ggplot(data = data.frame(K = seq(1, 100, by = 2), MSE = MSE_KNN), 
            aes(y = MSE, x = K)) +
       geom_line(alpha = 0.4, color = "darkgreen") +
       geom_point(alpha = 0.3, color = "darkgreen") +
       geom_hline(yintercept = MSE_ols, color = "darkgray", 
          alpha = 0.8, linetype = "dashed") +
       annotate("text", x = 90, y = 4e+13, label = "KNN", color = "darkgreen", alpha = 0.5) +
       annotate("text", x = 90, y = 3.7e+13, label = "OLS", color = "darkgray") +
       theme_bw()

     

     As observed, regardless of the choice of \(K\) in our KNN regression, the MSE obtained from KNN is consistently much smaller than that from OLS.

4. Curse of Dimensionality: When dealing with high-dimensional data (i.e., when you have a lot of predictors), KNN regression may suffer from the “curse of dimensionality.” In such cases, the performance of KNN deteriorates significantly as the dimensionality of the data increases. In the case of KNN regression, as the number of dimensions increases, the distance between data points becomes less meaningful. This is because in high-dimensional spaces, the concept of distance becomes less discriminating, as most data points are located far apart from each other. Consequently, KNN may struggle to find relevant neighbors, leading to poor predictive performance and increased computational complexity. On the other had, linear regression, with its fewer parameters, is less affected by this issue.

   

   Even in problems in which the dimension is small, we might prefer linear regression to KNN from an interpretability standpoint. If the test MSE of KNN is only slightly lower than that of linear regression, we might be willing to forego a little bit of prediction accuracy for the sake of a simple model that can be described in terms of just a few coefficients, and for which p-values are available.

2. Simple Linear Regression

Below is a population regression line equation, assuming that we are the god and know the data-generating process of the relationship between \(X\) (independent/explanatory variable) and \(Y\) (dependent/response/outcome variable).

However, since we usually don’t have the population data but only have access to a sample of it, we can only use the sample dataset to estimate the population regression line. In the below figure, \(b_0\) is the estimate of the regression intercept \(\beta_0\), while \(b_1\) is the estimate of the regression coefficient/slope \(\beta_1\).

The Ordinary Least Squares (OLS) approach aims to fit a line that minimizes the sum of the squared residuals, which means that the goal is to find a pair of \(b_0\) and \(b_1\) that can minimize the difference between our observed \(Y_i\) and estimated \(\hat{Y_i}\).

\[S(A, B) = min \Sigma E_i = min \Sigma (Y_i-\hat{Y_i})^2=min \Sigma (Y_i-(A+BX_i))^2\]

\[S(A, B) = \Sigma (Y_i-A-BX_i)^2\]

To find a pair of \(b_0\) and \(b_1\) that can minimize the difference between our observed \(Y_i\) and estimated \(\hat{Y_i}\) is an optimization problem. All we have to do is to take the partial derivatives of the above equation with respect to \(b_0\) and \(b_1\), respectively, set them equal to zero, and then solve it.

\[\frac{\partial S(A, B)}{\partial A}=\Sigma (-1)(2)(Y_i-A-BX_i)=0\]

\[\frac{\partial S(A, B)}{\partial B}=\Sigma (-X_i)(2)(Y_i-A-BX_i)=0\]

By solving the above two equations, you will get:

\[A = \bar{Y}-B\bar{X}\]

\[B = \frac{\Sigma (X_i-\bar{X})(Y_i-\bar{Y})}{\Sigma (X_i-\bar{X})^2}\]

a. Gauss-Markov Theorem

Gauss Markov theorem states that “if the errors are independently distributed with zero expectation and constant variance, then the least-squares estimator b is the most efficient linear unbiased estimator of \(\beta\). That is, of all unbiased estimators that are linear functions of the observations, the least-squares estimator has the smallest sampling variance and, hence, the smallest mean-squared error. For this reason, the least-squares estimator is sometimes termed BLUE, an acronym for best linear unbiased estimator” (Fox 2016: 212).

In other words, under certain conditions, the ordinary least squares (OLS) estimates are thought to be the best linear unbiased estimator (BLUE). Here, I review a set of assumptions (known as Gauss-Markov assumptions) that leads to a BLUE estimator.

1. Linearity: the expected (mean) value of the disturbance term is 0 (i.e., \(E(u_i)=0\)).

   Fox (2016) explains that a violation of this assumption implies “that the model fails to capture the systematic pattern of relationship between the response and explanatory variables” (307). This assumption could be broken in two ways: one, the mean of error could be a constant across all values of \(X_i\), two, it could vary. In the former case, the intercept term will be consistently off by the error amount. A uniform measurement error across all observations would be a cause of such an error (Berry 1993: 42). A common cause of the latter case will be omitted variable bias. Omitted variable bias deserves particular attention because it can cause the OLS estimate to be both biased and inconsistent. Another common cause is the wrong functional form (e.g., specifying a linear relation instead of a quadratic relation). Thus, when nonlinearity is detected, we should explore not only variable transformations and different model specifications (Fox, 2016: 317), but also the possibility of omitted variables. In short, before running OLS, reseachers should make sure that the inclusion of relevant variables in the model and the setting of the function are resonable and are based on the theoretical understanding.

2. Nonstochastic regressors: \(X\) values are independent of the error term (i.e., \(cov(X_i, u_i)=0\)).

   The violation of this assumption suggests that there is an endogeneity problem between the explanatory and outcome variables, which can arise from measurement error on \(X\), omitted confounders, or reciprocal causation between \(X\) and \(Y\). Again, if the main goal of the research is to provide an unbiased and consistent estimator to explain how a social phenomenon works in a real-life setting, to provide a convincing explanation requires us to make sure that the effect of our explanatory variable is actually exogenous.Fortunately, scholars have developed many methods to address this issue. For example, researchers could employ an instrumental variable or matching to improve the unbiasedness of the estimated effect.

3. Homoskedasticity: constant error variance across values of \(X_i\) (i.e., \(Var(\varepsilon)=\sigma^2\)).

   One instance where heteroskedasticity is often violated is cross-sectional studies (Berry, 1993: 73). For example, because more developed countries can have better data accuracy, measurement error could be correlated with the level of development. Thus, in a study that includes the level of development as one of the independent variables, the variance of error term may not be constant (Berry, 1993: 73). In such circumstances, Fox (2016) suggests that we can run weighted-least-squares (WLS) regression to overcome such probelm (304). Other solutions include transformation of the dependent variable and correction of coefficient standard errors (e.g. robust standard errors and clustered standard errors) for heteroskedasticity (Fox, 2016: 307).

4. Independence: the observations are sampled independently. no autocorrelation between disturbances (i.e., \(cov(u_i, u_j)=0\) for \(i \neq j\)).

   The independence assumption is important for time-series analysis, where the assumption of spherical errors is often violated. To overcome this, Beck & Katz (1995) propose panel-correlated standard errors to when analyzing time-series data.

5. Normality: errors are distributed normally (i.e.,\(\varepsilon \sim N(0, \sigma^2)\)).

Assumptions of 1 and 2 are related to the bias of the estimator, whereas assumptions of 3 and 4 are related to the efficiency of the estimator: when the assumption of 1 or 2 is violated, the OLS estimator is no longer unbiased; when the assumption of 3 or 4 is violated, the OLS estimator is no longer efficient. With these weak set of assumptions (1-4), according to the Gauss-Markov theorem, the OLS estimator is BLUE.

(In statistics, efficiency is a measure of the quality of an estimator, which can be characterized by having a smaller possible variance.)

With the fifth assumption, normality, the OLS estimator is the most efficient among all the unbiased estimators, not just linear estimators. Adding this fifth assumption makes the Gauss-Markov theorem a strong set. This last assumption is important when we aim to make a causal inference using the OLS estimator. With the normal distribution of the errors assumption, we can also infer that the OLS estimators also have normal sampling distributions. As a result, it allows us to apply statistical tests such as t-tests even for small sample size studies.

One important caveat is that when there is perfect multicollinearity (i.e. when two or more independent variables are perfectly correlated), we cannot even get an estimate using the OLS method. However, even when there is high multicollinearity, the OLS estimate is still BLUE (Berry 1993: 27). In such cases, standard errors will be very high, making the estimates fluctuate considerably from sample to sample (Berry 1993: 27).

The below table presents the consequences of the violation of Gauss-Markov assumptions and corresponded suggested solutions.

Assumption	Violation	Solution
Linearity	Biased/inconsistent estimates	Transformations, polynomials, different model
Nonstochastic regressors	Biased/inconsistent estimates	Instrumental variables
Homoskedasticity	Biased standard errors	Robust standard errors
Independence	Biased standard errors	Fixed effects
Perfect collinearity	cannot run OLS	Omit one collinear term

b. Properties of Least Squares Estimator

The sample least squares coefficients are unbiased estimators of the population regression coefficients.

Here are some reviews before demonstrating the unbiasedness of the least-squares estimators \(A\) and \(B\) for \(\alpha\) and \(\beta\):

• The assumed model: \(Y_i=\alpha+\beta X_i+\varepsilon_i\)
• Linearity assumption: \(E(\varepsilon_i)=0\)
• Normality assumption: \(\varepsilon_i\sim N(0, \sigma_\varepsilon^2)\)
• Independence assumption: \(Cov(\varepsilon_i, \varepsilon_j)=0, for i\neq j\)
• The expectation of linear combination: \(E(a+bY)=a+bE(Y)\)
• The variance of linear combination:
  • \(Var(a+bY)=b^2Var(Y)\)
  • \(Var(X+Y)=Var(X)+Var(Y)+2Cov(X, Y)\)
• For the below demonstration, the \(X_i\) are assumed to be fixed, not random.
  • \(B=\frac{\Sigma(X_i-\bar{X})(Y_i-\bar{Y})}{\Sigma(X_i-\bar{X})^2}=\frac{\Sigma(X_i-\bar{X})Y_i-\Sigma(X_i-\bar{X})\bar{Y}}{\Sigma(X_i-\bar{X})^2}=\frac{\Sigma(X_i-\bar{X})Y_i-\bar{Y}\Sigma(X_i-\bar{X})}{\Sigma(X_i-\bar{X})^2}=\frac{\Sigma(X_i-\bar{X})Y_i}{\Sigma(X_i-\bar{X})^2}\)
  • \(\Sigma(X_i-\bar{X})^2=\Sigma(X_i-\bar{X})(X_i-\bar{X})=\Sigma X_i(X_i-\bar{X})-\bar{X}\Sigma(X_i-\bar{X})=\Sigma X_i(X_i-\bar{X})\)

1. Demonstrate the unbiasedness of the least-squares estimators \(B\) for \(\beta\) in simple regression (i.e., \(E(B) = \beta\)).

   \(E(B)=\frac{\Sigma(X_i-\bar{X})(Y_i-\bar{Y})}{\Sigma(X_i-\bar{X})^2}\)

   \(=\frac{\Sigma(X_i-\bar{X})Y_i}{\Sigma(X_i-\bar{X})^2}\)

   \(=\frac{1}{\Sigma(X_i-\bar{X})^2}E(\Sigma(X_i-\bar{X})Y_i)\)

   \(=\frac{1}{\Sigma(X_i-\bar{X})^2}\Sigma(X_i-\bar{X})E(Y_i)\)

   \(=\frac{1}{\Sigma(X_i-\bar{X})^2}\Sigma(X_i-\bar{X})E(\alpha + \beta X_i+\varepsilon_i)\)

   \(=\frac{1}{\Sigma(X_i-\bar{X})^2}\Sigma(X_i-\bar{X})(\alpha +\beta X_i+E(\varepsilon_i))\)

   \(=\frac{1}{\Sigma(X_i-\bar{X})^2}\Sigma(X_i-\bar{X})(\alpha+\beta X_i)\)

   \(=\frac{1}{\Sigma(X_i-\bar{X})^2}(\Sigma(X_i-\bar{X})\alpha+\Sigma(X_i-\bar{X})\beta X_i)\)

   \(=\frac{\beta}{\Sigma(X_i-\bar{X})^2}(\Sigma(X_i-\bar{X})X_i)\)

   \(=\beta\)

2. Demonstrate the unbiasedness of the least-squares estimators \(A\) for \(\alpha\) in simple regression (i.e., \(E(A) = \alpha\)).

   \(E(A)=E(\bar{Y}-B\bar{X})\)

   \(=E(\frac{1}{n}\Sigma Y_i-\frac{\Sigma(X_i-\bar{X})(Y_i-\bar{Y})}{\Sigma(X_i-\bar{X})^2}\bar{X})\)

   \(=E(\frac{1}{n}\Sigma Y_i-\frac{\Sigma(X_i-\bar{X})Y_i}{\Sigma(X_i-\bar{X})^2}\bar{X})\)

   \(=E(\Sigma(\frac{1}{n}-\frac{\bar{X}(X_i-\bar{X})}{\Sigma(X_i-\bar{X})^2})Y_i)\)

   \(=E(\Sigma(\frac{1}{n}-\frac{\bar{X}(X_i-\bar{X})}{\Sigma(X_i-\bar{X})^2})(\alpha +\beta X_i +\varepsilon_i))\)

   \(=E(\Sigma(\frac{1}{n}-\frac{\bar{X}(X_i-\bar{X})}{\Sigma(X_i-\bar{X})^2}))E(\alpha +\beta X_i +\varepsilon_i)\)

   \(=E(\Sigma(\frac{1}{n}-\frac{\bar{X}(X_i-\bar{X})}{\Sigma(X_i-\bar{X})^2}))(\alpha+\beta X_i+E(\varepsilon_i))\)

   \(=\Sigma(\frac{1}{n}-\frac{\bar{X}(X_i-\bar{X})}{\Sigma(X_i-\bar{X})^2})(\alpha+\beta X_i)\)

   \(=\Sigma \frac{\alpha}{n}+\Sigma \beta \frac{X_i}{n}-\alpha \frac{\bar{X}\Sigma(X_i-\bar{X})}{\Sigma(X_i-\bar{X})^2}-\beta \frac{\bar{X}\Sigma(X_i-\bar{X})X_i}{\Sigma(X_i-\bar{X})^2}\)

   \(=\alpha + \beta \bar{X}-\alpha*0-\beta \bar{X}\)

   \(=\alpha\)

3. Derive the sampling variances of \(A\) in a simple regression.

   \(Var(B)=Var(\frac{\Sigma(X_i-\bar{X})Y_i}{\Sigma(X_i-\bar{X})^2})\)

   \(=\frac{1}{(\Sigma(X_i-\bar{X})^2)^2}Var(\Sigma(X_i-\bar{X})Y_i)\)

   \(=\frac{1}{(\Sigma(X_i-\bar{X})^2)^2}Var(\Sigma(X_i-\bar{X})(\alpha + \beta X_i+\varepsilon_i))\)

   \(=\frac{1}{(\Sigma(X_i-\bar{X})^2)^2}Var(\Sigma(X_i-\bar{X}(\alpha + \beta X_i)+\Sigma(X_i-\bar{X})\varepsilon_i)\)

   \(=\frac{1}{(\Sigma(X_i-\bar{X})^2)^2}Var(\Sigma(X_i-\bar{X})\varepsilon_i)=\frac{1}{(\Sigma(X_i-\bar{X})^2)^2}\Sigma Var((X_i-\bar{X})\varepsilon_i)\)

   \(=\frac{1}{(\Sigma(X_i-\bar{X})^2)^2}(\Sigma (X_i-\bar{X})^2 Var(\varepsilon_i))\)

   \(=\frac{1}{(\Sigma(X_i-\bar{X})^2)^2}(\Sigma (X_i-\bar{X})^2 \sigma_\varepsilon^2)\)

   \(=\frac{\sigma_\varepsilon^2}{\Sigma(X_i-\bar{X})^2}\)

4. Derive the sampling variances of \(B\) in a simple regression.

   \(Var(A)=Var(\bar{Y}-B\bar{X})\)

   \(=Var(\bar{Y})+Var(B\bar{X})-2Cov(\bar{Y},B\bar{X})\)

   \(=Var(\bar{Y})+Var(B\bar{X})\)

   \(=Var(\frac{\Sigma Y_i}{n})+\bar{X}^2Var(B)\)

   \(=\frac{1}{n^2}Var(\Sigma Y_i)+\bar{X}^2Var(B)\)

   \(=\frac{1}{n^2}\Sigma Var(Y_i)+\bar{X}^2Var(B)\)

   \(=\frac{1}{n^2}\Sigma Var(\alpha+\beta X_i+\varepsilon_i)+\bar{X}^2Var(B)\)

   \(=\frac{1}{n^2}\Sigma Var(\varepsilon_i)+\bar{X}^2Var(B)\)

   \(=\frac{1}{n^2}n\sigma_\varepsilon^2+\bar{X}(\frac{\sigma_\varepsilon^2}{\Sigma(X_i-\bar{X})^2})\)

   \(=\sigma_\varepsilon^2(\frac{1}{n}+\frac{\bar{X}^2}{\Sigma(X_i-\bar{X})^2})\)

   \(=\sigma_\varepsilon^2(\frac{\Sigma(X_i-\bar{X})^2 +n\bar{X}^2}{n\Sigma(X_i-\bar{X})^2})\)

   \(=\sigma_\varepsilon^2(\frac{\Sigma X_i^2-2\bar{X}\Sigma X_i+\Sigma \bar{X}^2+n\bar{X}^2}{n\Sigma(X_i-\bar{X})^2})\)

   \(=\sigma_\varepsilon^2(\frac{\Sigma X_i^2-2\bar{X}\Sigma X_i+2n \bar{X}^2}{n\Sigma(X_i-\bar{X})^2})\)

   \(=\sigma_\varepsilon^2(\frac{\Sigma X_i^2-2\frac{\Sigma X_i}{n}\Sigma X_i+2n \bar{X}^2}{n\Sigma(X_i-\bar{X})^2})\)

   \(=\sigma_\varepsilon^2(\frac{\Sigma X_i^2-\frac{2\Sigma X_i^2}{n}+2n\frac{\Sigma X_i^2}{n^2}}{n\Sigma(X_i-\bar{X})^2})\)

   \(=\frac{\sigma_\varepsilon^2 \Sigma X_i^2}{n\Sigma(X_i-\bar{X})^2}\)

c. Confidence Interval & Hypothesis Testing

As mentioned earlier, the fifth assumption, normality, allows us to apply statistical inference tests to assess the \(\beta\)s we estimate. Recall from the lecture that the residual standard error \(S_E\) (how closely the regression line we estimate fits the scatter of our data points) is:

\[\hat{\sigma } = SE(E_i) = \sqrt{\frac{\Sigma (E_i)^2}{n-2}} = \sqrt{\frac{\Sigma (Y_i-\hat{Y_i})^2}{n-2}}\]

And the standard error of the sample intercept (\(A\)) and slope (\(B\)) are:

\[SE(A) = \hat{\sigma }\sqrt{\frac{\Sigma X_i^2}{n\Sigma (X_i-\bar{X})^2}}\]

\[SE(B) = \frac{\hat{\sigma }}{\sqrt{\Sigma (X_i-\bar{X})^2}}\]

With standard errors, we can construct a \(100(1-\alpha)%\) confidence interval for our slope and perform hypothesis testing.

• Confidence Interval

  \(\beta=B+t_{\frac{\alpha}{2}}SE(B)\)

• Hypothesis Test

  Two-tailed test (by defualt, lm()function conducts two-tailedtest)

  \(H_0: \beta = 0\)

  \(H_1: \beta \neq 0\)

  \(t=\frac{b_1-\beta_1}{SE(b_1)}=\frac{b_1-0}{SE(b_1)}\)

  Compare this t-statistic with \(t_{\frac{\alpha}{2}, df=n-2}\) to see if it is larger or smaller than the critical value or not. If it is larger or smaller than the critical value, we can reject \(H_0\) in favor of \(H_1\). In addition to calculating t-statistics, we can also calculate p-value based the t-statistics and the critical values to perform hypothesis testing (p-value \(= 2*Pr(t \geq |t_c|)\) vs \(\alpha\)).

  One-tailed test (if your hypotheses are directional, you can also conduct one-tailed test)

  \(H_0: \beta = 0\)

  \(H_1: \beta < or > 0\)

  \(t=\frac{B-\beta_1}{SE(B)}=\frac{B-0}{SE(B)}\)

  When performing a one-tailed test, compare the t-statistic with \(t_{\alpha, df=n-2}\) without dividing \(\alpha\) by 2. P-value in one-tailed test is calculated as:

  p-value \(= Pr(t \geq or \leq t_c)\) vs \(\alpha\)

  Let’s run a simple example of it using the prestige data from Problem Set 1. The prestige data consists of 10 observations with 2 variables. The description of the variables are as follows:

  • prestige:The average prestige rating for the occupation.

  • education: The average number of years of education for occupational incumbents.

  # load data
  df <- read.csv("/Users/yu-shiuanhuang/Desktop/method-sequence/data/ps1_prestige.csv") 
  df

  ##    prestige education
  ## 1        82        86
  ## 2        83        76
  ## 3        90        92
  ## 4        76        90
  ## 5        90        86
  ## 6        87        84
  ## 7        93        93
  ## 8        90       100
  ## 9        52        87
  ## 10       88        86

  Let’s first calculate all the statistics we’ll need to perform hypothesis testing.

  1. Slope: \(B\)

     \(B=\frac{\Sigma(X_i-\bar{X})(Y_i-\bar{Y})}{\Sigma(X_i-\bar{X})^2}\)

  # calculate the slope b1
  B <- sum((df$education-mean(df$education))*(df$prestige-mean(df$prestige)))/
  sum((df$education-mean(df$education))^2)
  B

  ## [1] 0.3895028

  2. Intercept: \(A\)

     \(A=\bar{Y}-B\bar{X}\)

  # calculate the intercept b0
  A <- mean(df$prestige)-B*mean(df$education)
  A

  ## [1] 48.82376

  3. Residual Standard Error: \(\hat{\sigma} = SE(E_i)\)

     \(\hat{\sigma } = SE(E_i) = \sqrt{\frac{\Sigma (E_i)^2}{n-2}} = \sqrt{\frac{\Sigma (Y_i-\hat{Y_i})^2}{n-2}}\)

  # generate predicted Y
  df$pred <- A + B*df$education
  
  #calculate the residual standard error 
  se <- sqrt(sum((df$prestige-df$pred)^2)/(nrow(df)-2))
  se

  ## [1] 12.46986

  4. Standard Error of the Estimated Intercept and Slope: \(SE(A)\) & \(SE(B)\)

     \(SE(A) = \hat{\sigma }\sqrt{\frac{\Sigma X_i^2}{n\Sigma (X_i-\bar{X})^2}}\)

     \(SE(B) = \frac{\hat{\sigma }}{\sqrt{\Sigma (X_i-\bar{X})^2}}\)

  # calculate the standard error of the estimated intercept b0
  se_A <- se*sqrt(sum(df$education^2)/(nrow(df)*sum((df$education-mean(df$education))^2)))
  se_A

  ## [1] 57.80998

  # calculate the standard error of the estimated slope b1
  se_B <- se/sqrt(sum((df$education-mean(df$education))^2))
  se_B

  ## [1] 0.6554015

  Now, let’s perform two-tailed hypothesis testing of the slope.

  \(H_0: \beta = 0\)

  \(H_1: \beta \neq 0\)

  \(t=\frac{B-\beta_1}{SE(B)}=\frac{B-0}{SE(B)}=\frac{0.3895028-0}{0.6554015}=0.5942964\)

  # calculate the t-statistic
  t <- B/se_B
  t

  ## [1] 0.5942964

  What are the t critical values in both tails? Let’s say the level of significance, \(\alpha\), is 0.05 in this case.

  \(|t_{\frac{\alpha}{2}, df=n-2}|=|t_{\frac{0.05}{2}, df=10-2}|=2.306004\)

  # find t critical value
  qt(p = 0.025, df = 8, lower.tail = TRUE)

  ## [1] -2.306004

  qt(p = 0.025, df = 8, lower.tail = FALSE)

  ## [1] 2.306004

  Now, let’s compare our t-statistic of \(b_1\) with the t critical value.

  \(t=0.5942964 < |t_{\frac{0.05}{2}, df=8}|=2.306004\)

  We can also calculate the p-value.

  p-value \(= 2*Pr(t \geq 0.5942964)=0.5687352 > 0.05\)

  # find p-value
  2*pt(q = 0.5942964, df = 8, lower.tail = FALSE)

  ## [1] 0.5687352

  These all suggest that we cannot reject \(H_0\), which means that the positive effect of education on one’s prestige is not statistically significant.

  Let’s also try perform one-tailed hypothesis testing for this case.

  \(H_0: \beta = 0\)

  \(H_1: \beta > 0\)

  \(t=\frac{B-\beta_1}{SE(b_1)}=\frac{B-0}{SE(B)}=\frac{0.3895028-0}{0.6554015}=0.5942964\)

  What is the t critical values in the right tail?

  \(t_{\alpha, df=n-2}=t_{0.05, df=10-2}=1.859548\)

  # find t critical value
  qt(p = 0.05, df = 8, lower.tail = FALSE)

  ## [1] 1.859548

  Now again, let’s compare our t-statistic of \(B\) with the t critical value.

  \(t=0.5942964 < t_{0.05, df=8}=1.859548\)

  And the p-value in the one-tailed test is:

  p-value \(=Pr(t \geq 0.5942964)= 0.2843676 > 0.05\)

  # find p-value
  pt(q = 0.5942964, df = 8, lower.tail = FALSE)

  ## [1] 0.2843676

  In R, lm() does all the above calculation for us!

  # simple regression
  fit <- lm(prestige ~ education, data = df)
  summary(fit)

  ## 
  ## Call:
  ## lm(formula = prestige ~ education, data = df)
  ## 
  ## Residuals:
  ##      Min       1Q   Median       3Q      Max 
  ## -30.7105   0.3157   4.9580   5.6238   7.9525 
  ## 
  ## Coefficients:
  ##             Estimate Std. Error t value Pr(>|t|)
  ## (Intercept)  48.8238    57.8100   0.845    0.423
  ## education     0.3895     0.6554   0.594    0.569
  ## 
  ## Residual standard error: 12.47 on 8 degrees of freedom
  ## Multiple R-squared:  0.04228,    Adjusted R-squared:  -0.07743 
  ## F-statistic: 0.3532 on 1 and 8 DF,  p-value: 0.5687

d. Evaluating Least Squares Fit

We understand that there will be variability in the variable prestige across different occupations, with some having high prestige and others having low prestige. Some of this variability can be attributed to factors not included as variables in our regression model, such as economic conditions and other unmeasured variables. Additionally, differences in the variable education may explain some of the differences in prestige across occupations. For instance, occupations with a higher average number of years of education among incumbents tend to have higher prestige. However, it’s important to quantify how much of the variability in prestige can be explained by the variables included in our regression model, and how much remains unexplained by other factors that we have not accounted for.

1. Explained Variability: Regression Sum of Squares (RegSS)

   The Regression Sum of Squares (RegSS) (sometimes referred to as the Sum of Squares Explained, SSE) is a measure of the variability in the outcome variable that is explained by the explanatory variables, i.e. the x-variables in your regression. It is given by the following sum:

   \(RegSS = \sum_{i=1}^{n}(\hat{y_i}-\bar{y})^2\)

2. Residual or Unexplained Variability: Residual Sum of Squares (RSS)

   The Residual Sum of Squares (RSS) is a measure of the variability in the outcome variable that is not explained by your regression, and therefore is due to all the other factors that affect prestige besides education. It is given by the following sum:

   \(RSS = \sum_{i=1}^{n}(y_i-\hat{y_i})^2 = \sum_{i=1}^{n}e_i^2\)

3. Total Variability: Total Sum of Squares (TSS)

   You can show mathematically that \(RSS + RegSS\) is equal to the following expression, which is referred to as the Total Sum of Squares (TSS):

   \(TSS = \sum_{i=1}^{n}(y_i-\bar{y_i})^2 = \sum_{i=1}^{n}(y_i-\hat{y_i})^2 + \sum_{i=1}^{n}(\hat{y_i}-\bar{y})^2 = RSS + RegSS\)

4. Coefficient of Determination: R-Squared value

   The coefficient of determination, sometimes referred to as the R-Squared value, is a measure of what percentage of the variability in your outcome variable is explained by your explanatory variables. It is given by the expression,

   \(R^2 = \frac{RegSS}{TSS}\)

Discussion II

For today’s discussion, we will review concepts regarding how to find a set of \(\hat{\beta}\)s that can minimize the squared residual errors of our multivariate linear model using matrix algebra.

1. Multivariate Linear Regression
2. Find Coefficients using Matrix in R

1. Multivariate Linear Regression

The central difference between simple and multiple linear regression is that the slope coefficients for the explanatory variables in multiple regression are partial coefficients. That is, it represents the “effect” on the response variable of a one-unit increment in the corresponding explanatory variable, holding constant the value of the other explanatory variable.

\[Y_i=\hat{\beta_0}+\hat{\beta_1}X_{i1}+\hat{\beta_2}X_{i2}+...+\hat{\beta_k}X_{ik}+\varepsilon_i\]

Recall that the goal of the OLS approach is to minimize the discrepancy between our observed \(Y_i\) and estimated \(\hat{Y_i}\). To find a set of estimated \(\beta\)s that can minimize this discrepancy is an optimization problem. All we have to do is to take the partial derivatives of the above equation with respect to each coefficient, set them equal to zero, and then solve it. However, as shown in the lecture, adding more and more \(X\) variables to our linear regression model makes the calculation more and more tedious. In order to find the estimated \(\beta\)s more efficiently, we use matrix algebra to help us do the calculation.

To do this, let’s first think of our data set as a matrix. The below equations are how each observation in our data set generated.

\[Y_1=\hat{\beta_0}+\hat{\beta_1}X_{11}+\hat{\beta_2}X_{12}+...+\hat{\beta_k}X_{1k}+\varepsilon_1\]

\[Y_2=\hat{\beta_0}+\hat{\beta_1}X_{21}+\hat{\beta_2}X_{22}+...+\hat{\beta_k}X_{2k}+\varepsilon_2\]

.

.

.

\[Y_n=\hat{\beta_0}+\hat{\beta_1}X_{n1}+\hat{\beta_2}X_{n2}+...+\hat{\beta_k}X_{nk}+\varepsilon_n\]

Rewrite them into a matrix form.

The system of equations is summarized as the y vector equal to the product of the X matrix of data times the \(\hat{\beta}\) vector of parameters plus the vector of disturbances, \(\varepsilon\). That said, the fitted linear model we are trying to find is: \(y = X\hat{\beta}+\varepsilon\). In the lecture, Lauren showed how to apply what we did when finding a pair of intercept and slope in simple linear regression to solve the optimization problem(i.e., minimize RSS) to find the coefficient vector \(\hat{\beta}\) using matrix algebra.

\[ \begin{aligned} S(\hat{\beta}) &= \sum E_i^2 = \varepsilon\varepsilon' = (y- X\hat{\beta})'(y- X\hat{\beta})\\ &= y'y-y'X\hat{\beta}-\hat{\beta}'X'y+\hat{\beta}'X'X\hat{\beta} \\ &= y'y-(2y'X)\hat{\beta}+\hat{\beta}'(X'X)\hat{\beta} \end{aligned} \]

Instead of taking derivative with respect to each coefficient one at a time, we can take derivative directly with respect to the coefficient vector \(\hat{\beta}\).

\[ \begin{aligned} S(\hat{\beta}) &= \sum E_i^2 = \varepsilon\varepsilon' = (y- X\hat{\beta})(y- X\hat{\beta})'\\ &= y'y-y'X\hat{\beta}-\hat{\beta}'X'y+\hat{\beta}'X'X\hat{\beta} \\ &= y'y-(2y'X)\hat{\beta}+\hat{\beta}'(X'X)\hat{\beta} \end{aligned} \] \[ \begin{aligned} \partial S(\hat{\beta}) &= \frac{\partial}{\partial \hat{\beta}}(y'y-(2y'X)\hat{\beta}+\hat{\beta}'(X'X)\hat{\beta}) \\ 0 &=0-2Xy'+ 2\hat{\beta}X'X \\ 2\hat{\beta}X'X &=2Xy'\\ \hat{\beta} &= \frac{Xy'}{X'X} \\ \hat{\beta} &= (X'X)^{-1}Xy' \end{aligned} \] Therefore, the least squares estimator is:

\[\hat{\beta}=(X'X)^{-1}X'y\]

2. Find Coefficients using Matrix in R

Let’s use the Anscombe.txt data set as an example to find the estimated coefficients by using the above matrix equation.

• education = per capita education expenditures, dollars
• income = per capita income, dollars
• under18 = proportion under 18 years old, per 1000
• urban = proportion urban, per 1000

# load data
df2 <- read.delim("/Users/yu-shiuanhuang/Desktop/method-sequence/data/Anscombe.txt", sep = "")
head(df2)

##    education income under18 urban
## ME       189   2824   350.7   508
## NH       169   3259   345.9   564
## VT       230   3072   348.5   322
## MA       168   3835   335.3   846
## RI       180   3549   327.1   871
## CT       193   4256   341.0   774

# convert each variable to a vector
edu <- df2$education
intercept <- rep(1, nrow(df2)) # remember to crease this vector with 1s!
inc <- df2$income
un18 <- df2$under18
urb <- df2$urban
X <- cbind(intercept, inc, un18, urb)
class(X) # this is a matrix

## [1] "matrix" "array"

We now have a vector edu and a matrix \(X\) containing our intercepts and our three explanatory variables (income, under18, urban). Let’s estimate our best fit parameters \(\beta_0, \beta_1,\), \(\beta_2\), and \(\beta_3\) by applying the estimator we derived in the lecture: \(\hat{\beta} = (X'X)^{-1}X'y\).

beta.hat <- solve(t(X) %*% X) %*% t(X) %*% edu

t(beta.hat)

##      intercept        inc      un18        urb
## [1,] -286.8388 0.08065325 0.8173377 -0.1058062

Compare with the canned routine in R to confirm the matrix calculation was correct.

mod <- lm(edu ~ X-1) # remember to remove the first column in the X matrix bc that's the column with all the 1s we created just for the algebra calculation

beta.canned <- mod$coefficients
beta.canned

##    Xintercept          Xinc         Xun18          Xurb 
## -286.83876273    0.08065325    0.81733774   -0.10580623

How to come up with the predicted \(\hat{y}\)? \[\hat{y}=Xb\]

eduhat <- beta.hat[1] + df2$income*beta.hat[2] + 
    df2$under18*beta.hat[3] + df2$urban*beta.hat[4]
eduhat

##  [1] 173.8168 199.0526 211.7006 207.0077 174.5936 253.2396 223.9983 210.5846
##  [9] 179.8788 206.2476 213.3513 231.5918 233.2380 209.4856 211.0236 196.9737
## [17] 176.3859 188.1518 199.2828 195.3128 187.4341 250.0855 231.5108 252.0303
## [25] 182.3619 139.8510 169.8280 162.6433 176.5621 159.9772 156.6485 139.1353
## [33] 135.8968 148.7554 135.1561 174.0988 143.0524 175.0569 195.4585 171.6380
## [41] 205.2510 192.1739 197.6281 189.7953 190.1857 261.4657 212.7525 181.6204
## [49] 221.7760 355.7228 221.5298

How to come up with the residuals? \[e=y-\hat{y}=y-Xb\]

e <- edu - eduhat
e

##  [1]  15.1831901 -30.0526055  18.2993676 -39.0077435   5.4064124 -60.2396370
##  [7]  37.0016640   3.4153951  21.1211678 -34.2476465 -19.3513498 -42.5918115
## [13]  -0.2380450  -0.4855509  50.9763644  37.0263340   0.6140641 -11.1518234
## [19] -12.2828447 -47.3127738   8.5658819  -2.0854980  15.4891839  -6.0302784
## [25]  -2.3619147   9.1490037 -14.8279988 -13.6433457 -20.5621494  31.0227604
## [31] -16.6485112  -2.1352962 -23.8967859 -18.7553867  -1.1561315 -12.0987770
## [37]  -8.0523620 -20.0569429  42.5415211  -1.6380384  32.7489708  -0.1738631
## [43]  29.3718741  17.2047442  10.8143055 -36.4656906   2.2475110  51.3796304
## [49]  51.2240417  16.2771794  -9.5297657

Now calculate the standard errors of the vector \(\hat{\beta}\). Recall the error variance is unknown but we have available the unbiased estimator \(S^2_E = (e'e)/(n-k-1)\), based on the residuals of our model fit. We can extract these residuals, calculate the estimator, and then use it to calculate the variance of the least squares coefficients \(\hat{V(b)} = S^2_E (X'X)^{-1}\).

• As a quick refresher of concepts: the variance is a measure of a random variable’s “spread” or variation around its mean (a.k.a. its expected value), while the co-variance measures how correlated are the variations of two random variables with each other.

• The variance-covariance matrix is a square matrix (i.e. it has the same number of rows and columns). The elements of the matrix that lie along its main diagonal (i.e. the one that goes from top-left to bottom-right contain the variances while all other elements contain the co-variances). Thus, the variance-covariance matrix of the fitted coefficients of a regression model contains the variances of the fitted model’s coefficient estimates along its main diagonal, and it also contains the pair-wise co-variances between coefficient estimates in the non-diagonal elements.

• Why do we not only calculate the variance of regression coefficients but also the convariance of them?

  The covariance matrix can be very useful for certain model diagnostics. If two variables are highly correlated, one way to think about it is that the model is having trouble figuring out which variable is responsible for an effect (because they are so closely related). This can be helpful for a whole variety of cases, such as choosing subsets of covariates to use in a predictive model; if two variables are highly correlated, you may only want to use one of the two in your predictive model.

S2E <- as.numeric((t(e)%*%e)/(nrow(df2)-3-1))
S2E

## [1] 712.5394

v.beta.hat <- S2E*solve(t(X)%*%X)
v.beta.hat # the diagonal numbers are the variance of each coefficient

##              intercept           inc          un18           urb
## intercept 4214.5975090 -1.849526e-01 -9.7493677311 -0.1582902773
## inc         -0.1849526  8.646281e-05  0.0001387409 -0.0002162611
## un18        -9.7493677  1.387409e-04  0.0255327479  0.0002084764
## urb         -0.1582903 -2.162611e-04  0.0002084764  0.0011752674

sqrt(diag(v.beta.hat)) 

##    intercept          inc         un18          urb 
## 64.919931523  0.009298538  0.159789699  0.034282173

Compare with the canned routine in R to see if we got the same standard errors of coefficients.

	OLS
Variables	Estimates	S.E.	C.I. (95%)
Xintercept	-286.84 ***	64.92	-417.44 – -156.24
Xinc	0.08 ***	0.01	0.06 – 0.10
Xun18	0.82 ***	0.16	0.50 – 1.14
Xurb	-0.11 **	0.03	-0.17 – -0.04
Observations	51
R2 / R2 adjusted	0.984 / 0.982
p<0.05   ** p<0.01   *** p<0.001

With the standard errors of each coefficient, we can perform hypothesis testing and construct confidence intervals to determine the significance and reliability of the estimated effects (similar to what we did in the case of simple linear regression).

• What would happen when there is perfect multicollinearity in our data (i.e. when two or more independent variables are perfectly correlated)?

  # convert each variable to a vector
  edu <- df2$education
  intercept <- rep(1, nrow(df2)) # remember to crease this vector with 1s!
  inc <- df2$income
  un18 <- df2$under18
  urb <- df2$urban
  pmc <- 2*inc
  cor(pmc, inc) # pmc is perfectly correlated with inc

  ## [1] 1

  X <- cbind(intercept, inc, un18, urb, pmc)
  class(X) # this is a matrix

  ## [1] "matrix" "array"

  beta.hat <- solve(t(X) %*% X) %*% t(X) %*% edu

  ## Error in solve.default(t(X) %*% X): Lapack routine dgesv: system is exactly singular: U[5,5] = 0

  As you can see in the result of the above chunk, you cannot calculate the matrix algebra when there is perfect multicollinearity. However, when you run lm() in R, R will automatically delete the variable that is perfectly correlated with the other variable(s) in the dataset.

  mod2 <- lm(edu ~ X-1) # remember to remove the first column in the X matrix bc that's the column with all the 1s we created just for the algebra calculation

  	OLS
  Variables	Estimates	S.E.	C.I. (95%)
  Xintercept	-286.84 ***	64.92	-417.44 – -156.24
  Xinc	0.08 ***	0.01	0.06 – 0.10
  Xun18	0.82 ***	0.16	0.50 – 1.14
  Xurb	-0.11 **	0.03	-0.17 – -0.04
  Observations	51
  R2 / R2 adjusted	0.984 / 0.982
  p<0.05   ** p<0.01   *** p<0.001

  If there is high multicollinearity between variables but not perfectly correlated, the matrix algebra still works, and lm() will still spit out the estimates of the highly correlated variable. However, the standard errors of coefficients will be very high.

  # convert each variable to a vector
  edu <- df2$education
  intercept <- rep(1, nrow(df2)) # remember to crease this vector with 1s!
  inc <- df2$income
  un18 <- df2$under18
  urb <- df2$urban
  set.seed(1234)
  pmc <- 2*inc + rnorm(nrow(df2), 5, 1)
  cor(pmc, inc) # pmc is highly but not perfectly correlated with inc

  ## [1] 0.9999997

  X <- cbind(intercept, inc, un18, urb, pmc)
  class(X) # this is a matrix

  ## [1] "matrix" "array"

  beta.hat <- solve(t(X) %*% X) %*% t(X) %*% edu
  
  t(beta.hat)

  ##      intercept     inc      un18        urb       pmc
  ## [1,] -290.6956 -1.4651 0.8188972 -0.1038287 0.7726427

  eduhat <- beta.hat[1] + df2$income*beta.hat[2] + 
  df2$under18*beta.hat[3] + df2$urban*beta.hat[4] + pmc*beta.hat[5]
  e <- edu - eduhat
  S2E <- as.numeric((t(e)%*%e)/(nrow(df2)-4-1))
  v.beta.hat <- S2E*solve(t(X)%*%X)
  sqrt(diag(v.beta.hat))

  ##   intercept         inc        un18         urb         pmc 
  ## 69.35278085  9.01847688  0.16172179  0.03651251  4.50787061

  	OLS
  Variables	Estimates	S.E.	C.I. (95%)	Estimates	S.E.	C.I. (95%)
  Xintercept	-286.84 ***	64.92	-417.44 – -156.24	-290.70 ***	69.35	-430.30 – -151.10
  Xinc	0.08 ***	0.01	0.06 – 0.10	-1.47	9.02	-19.62 – 16.69
  Xun18	0.82 ***	0.16	0.50 – 1.14	0.82 ***	0.16	0.49 – 1.14
  Xurb	-0.11 **	0.03	-0.17 – -0.04	-0.10 **	0.04	-0.18 – -0.03
  Xpmc				0.77	4.51	-8.30 – 9.85
  Observations	51			51
  R2 / R2 adjusted	0.984 / 0.982			0.984 / 0.982
  p<0.05   ** p<0.01   *** p<0.001

Discussion III

For today’s discussion, we will review some diagnostic techniques to detect influential data points that could affect our model and ensure our model adheres to Gauss-Markov assumptions, addressing any violations as needed.

1. Unusual and Influential Data
2. Diagnostics for Errors and Linearity

1. Unusual and Influential Data

When dealing with observational or experimental data, it’s common to encounter outliers, which are data points that behave differently than the rest of the data for some reason. Outliers can greatly influence the analysis and interpretation of the data. In today’s discussion, we’ll examine what outliers are, how they can affect data analysis and interpretation, and use simple linear regression as an example to explore situations where individual data points significantly impact the model. Lastly, we will also

1.1 Leverage, discrepancy and influence

Although many researchers discuss outliers in a broad sense, it’s important to recognize that there are different types of unusual data points. A data point may be unusual in terms of its predictor behavior (the \(x\)-value in simple regression), its outcome (the \(y\)-value in simple regression), or both.

• leverage: A data point that differs significantly from the rest of the dataset in terms of its predictor behavior is known as a leverage point. In the case of simple linear regression, this means that the x-value of the data point is much higher or lower than the mean of the predictor variable.

• discrepancy: When a data point has an uncommon y-value for its corresponding x-value, it is said to have high discrepancy. Within the context of regression analysis, such an unusual point is called an outlier.

• influence: Having either a high leverage or a high discrepancy does not necessarily make a data point influential in a linear model. In fact, a single data point’s influence is determined by leverage × discrepancy. Therefore, having only high leverage or high discrepancy may not be enough to significantly alter the model parameters.

The figure below illustrates the three characteristics discussed above. In each of the three panels, the red line represents the line of best fit without the data point in question (marked by the triangle), while the blue line shows the line of best fit with it.

• In panel A, the data point with the triangle has a high leverage - its \(x\)-value is much higher than the rest.

• In panel B, it has a high discrepancy - it lies far away from the line of best fit.

  However, neither of the above characteristics alone exerts a lot of influence on our model parameters, as the red line does not diverge significantly from the blue line.

• In panel C, we see a data point that has both high leverage and high discrepancy, and as a result, it exerts high influence: the blue line is noticeably different from the red line.

Example Data

For today’s discussion, we will use a data set based on an example in Field, Miles, and Field (2012). The dataset includes eight samples, where the \(x\)-value represents the number of pubs within a borough (district) of London, and the \(y\)-value represents the number of deaths in that borough over a specific period of time. We are interested in investigating the relationship between the number of deaths and the number of pubs in each borough.

pubs <- data.frame(name = c(1, 2, 3, 4, 5, 6, 7, 8), 
                   pubs = c(10, 25, 40, 55, 70, 85, 100, 500), 
                   deaths = c(1043.822, 2086.934, 2951.086, 3992.459, 5088.003, 
                              6095.645, 6923.497, 10000.000))
pubs

##   name pubs    deaths
## 1    1   10  1043.822
## 2    2   25  2086.934
## 3    3   40  2951.086
## 4    4   55  3992.459
## 5    5   70  5088.003
## 6    6   85  6095.645
## 7    7  100  6923.497
## 8    8  500 10000.000

Let’s plot the line of best fit of this data:

library(tidyverse)
ggplot(data = pubs, aes(x = pubs, y = deaths, label = name)) +
    geom_point(alpha = 0.8) +
    geom_text(hjust = 2, vjust = 0) +
    geom_smooth(method = lm, se = FALSE, colour = "darkorange") +
    xlab("Number of Pubs") + 
    ylab("Deaths") +
    theme_bw()

Both the raw data and the linear regression plot indicate that data point 8 is significantly different from the other samples. This particular data point not only has a much higher number of pubs compared to the others, but it also appears to have a distinct relationship between the number of pubs and the number of deaths, unlike the other boroughs.

Since we have some concerns about this specific data point, let’s determine whether it exhibits a) high leverage, b) high discrepancy, and c) high influence.

1.2 Assessing Leverage

Recall that leverage quantifies the degree to which a predictor value differs from the other predictor values. In the case of simple linear regression, we can compute the distance between each predictor point (\(X_i\)) and the mean predictor value (\(\bar{X}\)) to measure leverage.

A standardized version of this distance is called hat-value and denoted by \(h_i\) (for \(i = {1,..., n}\)):

\[h_i = \frac{1}{n}+\frac{(X_i-\bar{X})^2}{\sum_{j=1}^{n}(X_j-\bar{X})^2}\]

The average hat value (\(\bar{h}\)) is defined as \(\frac{k+1}{n}\), in which \(k\) is the number of predictors and \(n\) the number of observations. Values of \(h\) are bound between \(\frac{1}{n}\) and 1, with 1 denoting highest leverage (highest distance from mean).

Based on our sample data, it’s obvious that data point 8 has the highest leverage among all the points. Let’s apply the formula above to compute the leverage value (\(h_8\)) for this data point.

# number of cases
n = nrow(pubs)

# distance to mean of point 8
numerator = (pubs$pubs[8] - mean(pubs$pubs))^2

# distance to mean of all the other points
denominator = sum((pubs$pubs - mean(pubs$pubs))^2)

# putting it together
h_8 = 1/n + numerator/denominator

h_8

## [1] 0.969302

The resulting hat-value is 0.969302. This is quite high - in fact, it’s very close to 1, the highest possible value.

You don’t need to manually calculate all of the hat values because R offers a convenient hatvalues function that can be applied to any linear model. To take advantage of this feature, we will begin by fitting a simple linear model using lm and then extracting the hat value of each data point.

# fitting linear model
mod <- lm(deaths ~ pubs, data = pubs)

# getting hatvalues and plotting them
plot(hatvalues(mod))

Remember that leverage alone does not mean a point exerts high influence, but it certainly means it’s worth investigating. Hat values are open to interpretation, but a cut-off value that is common is twice or three times the average \(\bar{h}\), meaning anything above that value should be looked at closer.

Let’s add a cut-off line into the above hatvalue figure.

# the average of hat-values
mean_h <- (1+1)/8

# adding a cut-off line
plot(hatvalues(mod))
abline(h = 2*mean_h, lty = 2)
abline(h = 3*mean_h, lty = 2)

As illustrated in the figure above, we can see that the value of \(h_8\) is greater than twice and three times the mean hat values (\(\bar{h}\)). Based on this observation, we can safely conclude that the value of \(h_8\) is highly unusual. In this case, \(h_8\) is definitely unusual!

After having assessed leverage, let’s look at discrepancy. How unusual is the \(y\)-value given its \(x\)-value?

1.3 Assessing Discrepancy

As we learned earlier, data points that do not fit well to the linear regression line are referred to as outliers, or high-discrepancy points. Typically, we assess the fit of the regression line using residuals, which measure the distance between a predicted value and the actual value. While we already know that point 8 has high leverage, it is worth noting that the line of best fit is quite close to its predicted value (see the left panel). Additionally, when we examine the residuals, point 8 does not appear to be further from the regression line than any other points (see the right panel).

Instead, we can look at standardized or studentized residuals.

• Standardized residuals are formed by calculating:

  \[E_i'=\frac{E_i}{S_E\sqrt{1-h_i}}\]

  where \(S_E=\sqrt{\frac{\sum E_i^2}{n-k-1}}\)

  Let’s calculate the standardized residual for our example data by hand (the 8th data point is the one we’re interested in):

  # residual for data point in original model
  Ei <- as.numeric(residuals(mod)[8])
  
  # estimate of sigma (standard deviation) for residuals
  S_E <- summary(mod)$sigma
  
  # hatvalue for point 8
  hi <- as.numeric(hatvalues(mod)[8])
  
  # putting it together
  Eprime <- Ei/(S_E*sqrt(1-hi))
  Eprime

  ## [1] -2.447433

  You can also just let R do the maths by calling rstandard on the original model:

  rstandard(mod)

  ##           1           2           3           4           5           6 
  ## -1.44154418 -0.89769067 -0.48020238  0.04464274  0.59853447  1.09293034 
  ##           7           8 
  ##  1.47205817 -2.44743280

  Note that the last value is the same as the one we calculated by hand! Let’s compare the residuals and the studentized residuals in plots side by side:

  

  The standardized residuals provide a better measure of discrepancy and can reveal outliers that were not detected by the normal residual plot. In our case, by looking at the standardized residual plot below, we can confirm that point 8 has a high discrepancy compared to the other points.

• Studentized residuals are calculated by fitting a model without the case for which the residual is calculated, and then scaling the resulting residual (\(E_i\)) by an estimate of the standard deviation of the residuals (\(S_{E(-i)}\)) and the point’s hat value (\(h_i\)):

  \[E_i^*=\frac{E_i}{S_{E(-i)}\sqrt{1-h_i}}\]

  Let’s calculate the studentized residual for our example data by hand (the 8th data point is the one we’re interested in):

  # model excluding point 8
  mod.red <- lm(deaths ~ pubs, data = pubs[-8,]) 
  
  # residual for data point in original model
  Ei <- as.numeric(residuals(mod)[8])
  
  # estimate of sigma (standard deviation) for residuals
  S_E <- summary(mod.red)$sigma
  
  # hatvalue for point 8
  hi <- as.numeric(hatvalues(mod)[8])
  
  # putting it together
  Estar <- Ei/(S_E*sqrt(1-hi))
  Estar

  ## [1] -54.5284

  You can also just let R do the maths by calling rstudent on the original model:

  rstudent(mod)

  ##            1            2            3            4            5            6 
  ##  -1.62765333  -0.88075356  -0.44703731   0.04075983   0.56346508   1.11482846 
  ##            7            8 
  ##   1.68127219 -54.52839824

  Note that the last value is the same as the one we calculated by hand! Let’s compare the residuals and the studentized residuals in plots side by side:

  

  The studentized residuals reveal clearly that point 8 has a high discrepancy, while this was not possible to see from the normal residual plots. You can also perform Bonferroni adjustment to test the statistical significance of the \(E_8^*\) to check if we can confidently say that point 8 is an outlier. Observations are considered outliers if their Bonferroni p is less than .05.

  library(car)
  outlierTest(mod)

  ##   rstudent unadjusted p-value Bonferroni p
  ## 8 -54.5284         3.9231e-08   3.1385e-07

  According to the above test, we can say that point 8 can be considered an outlier (high discrepancy).

1.4 Assessing Influence

As mentioned earlier, having high leverage or discrepancy alone does not necessarily mean a data point will have a high influence on the linear regression model. However, since our point of interest (point 8) has both high leverage and discrepancy, it is likely to have high influence on the model.

We’ll use something called Cook’s D, which relies on the standardized residuals to assess influence. A similar measure is DFFITS, which instead is based on studentized residuals. Another measure we won’t go into are DFBETAS, which measures the influence on each individual parameter, instead of the overall model.

Observations that combine high leverage with large studentized residual exert substantial influence on regression coeffcients. Cook’s D statistic provides a summary index of influence.

\(D_i=\frac{E_i'^2}{k+1}\times \frac{h_i}{1-h_i}\)

where the first term, the standardized residual, is a measure of discrepancy and the second is a measure of leverage.The higher the leverage and residuals, the higher the Cook’s distance.

Several interpretations for Cook’s distance exist. There isn’t a universally accepted rule for cut off points.

• One interpretation is to investigate any point over \(\frac{4}{n}\), where \(n\) is the number of observations.
• Other scholars suggest that any “large” \(D_i\) should be investigated. How large is “too large”? The consensus seems to be that a \(D_i\) value of more than 1 indicates an influential value, but you may want to look at values above 0.5.

Let’s use influencePlot() in the car package to check Cook’s D statistics!

##      StudRes       Hat      CookD
## 1  -1.627653 0.1813863  0.2302244
## 7   1.681272 0.1256287  0.1556728
## 8 -54.528398 0.9693020 94.5671657

1.5 Added-variable Plot (AV Plot)

All the above techniques can also be applied to multivariate linear regression models. Additionally, we introduce another very useful influence graph for visually detecting outliers in multivariate linear models, known as an added-variable plot or partial-regression plot.

A limitation of using the typical residual plots in a multivariate linear regression model, such as plotting residuals against the values of a predictor variable, is that they may not fully illustrate the “additional contribution” of a predictor variable when considering other variables already in the model. Added-variable plots offer a more refined visualization, providing graphic insights into the marginal importance of a predictor variable, taking into account the presence of other variables in the model. The following materials are sourced from Fox’s slides and lecture.

• Let \(y_i^{(1)}\) represent the residuals from the least-squares regression of \(y\) on all of the \(x\)s with the exception \(x_1\), that is, the residuals from the fitted model:

  \[y_i = b_0^{(1)} + b_2^{(1)}x_{i2} + ...+ b_k^{(1)}x_{ik} + y_i^{(1)}\]

• Likewise, \(x_i^{(1)}\) are residuals from the least-squares regression of \(x_1\) on the other \(x\)s:

  \[x_{1i} = c_0^{(1)} + c_2^{(1)}x_{i2} + ...+ c_k^{(1)}x_{ik} + x_i^{(1)}\]

• The notation emphasizes the interpretation of the residuals \(y_i^{(1)}\) and \(x_i^{(1)}\) as the parts of \(y\) and \(x_1\) that remain when the linear dependence of these variables on \(x_2, ..., x_k\) is removed.

The AV plot for \(x_1\) is then the scatterplot of \(y_i^{(1)}\) versus \(x_i^{(1)}\), and we repeat the procedure for each \(x_j\), \(j =0, 1, ..., k\) (where \(x_0 = 1\)). In effect, the \((k+1)\)-dimensional scatterplot for \(y\) and \(x_1, ..., x_k\) is reduced to a sequence of \(k+1\) 2D AV plots. By doing so, the AV plots visualize leverage and influence on each of the regression coefficients.

The added-variable plot for \(x_1\) has the following very interesting properties:

• The slope of the least-squares simple-regression line of of \(y_i^{(1)}\) on \(x_i^{(1)}\) is the same as the least-squares slope \(b_1\) for \(x_1\) in the full multiple regression.

• The residuals from this simple regression are the same as the residuals \(e_i\) from the full regression.

• Consequently, the standard deviation of the residuals in the added-variable plot is the same as the multiple regression.

• The standard error of \(b_1\) from this simple regression is also the same as the standard error from the full regression

• Because the \(x_i^{(1)}\) are residuals, they are less variable than \(x_1\) if \(x_1\) is correlated with the other \(x\)s. The added-variable plot therefore show how collinearity can degrade the precision of the estimation by decreasing the conditional variation of an \(x\).

NOTE: You may find these concepts very familiar as in POL 212, we have actually covered the Frisch-Waugh-Lovell Theorem (you can review the materials of Discussion IV in POL 212 Discussion), in which we showed how to decompose a regression of \(Y\) on a set of variables \(X\) into two pieces. The coefficient of one of the variables in multiple linear regression model can be obtained by netting off effect of other variables in the regression model from both dependent and independent variables. The FWL Theorem shows that it is important to partial one set of covariates out of the other in order to obtain the ceteris paribus effect. For example, in order to get the effect of \(X1\) on \(Y\), \(X2\) needs to be held constant, which involves removing the portion of \(X1\) that is correlated with \(X2\). Likewise, when identifying outliers among explanatory variables in a multivariate linear regression model, it’s crucial to first hold all other variables constant.

Example

For this example, we will use the Duncan data frame, which has 45 rows and 4 columns. This is a dataset on the prestige and other characteristics of 45 U.S. occupations in 1950.

library(car)
library(carData)

knitr::kable(head(Duncan)) # from carData package

	type	income	education	prestige
accountant	prof	62	86	82
pilot	prof	72	76	83
architect	prof	75	92	90
author	prof	55	90	76
chemist	prof	64	86	90
minister	prof	21	84	87

These are the description of each variable:

• type: Type of occupation. A factor with the following levels: prof, professional and managerial; wc, white-collar; bc, blue-collar.

• income: Percentage of occupational incumbents in the 1950 US Census who earned $3,500 or more per year (about $36,000 in 2017 US dollars).

• education: Percentage of occupational incumbents in 1950 who were high school graduates (which, were we cynical, we would say is roughly equivalent to a PhD in 2017).

• prestige: Percentage of respondents in a social survey who rated the occupation as “good” or better in prestige.

Let’s begin by regressing prestige on all explanatory variables income and education, then plot the added-variable plots to visually inspect for any potential outliers.

m.duncan <- lm(prestige ~ income + education, data = Duncan)

	OLS
Variables	Estimates	S.E.	C.I. (95%)
(Intercept)	-6.06	4.27	-14.69 – 2.56
income	0.60 ***	0.12	0.36 – 0.84
education	0.55 ***	0.10	0.35 – 0.74
Observations	45
R2 / R2 adjusted	0.828 / 0.820
p<0.05   ** p<0.01   *** p<0.001

avPlots(m.duncan)

Note that the slope of the blue regression line in both plots corresponds to the estimated coefficient of income and education in our multivariate regression model.

In the AV plots, we can observe the impact (or lack thereof) of these unusual cases on the regression coefficients. Let’s first examine the left-hand side of the figure. We notice that the case of the minister pulls up on the left side, while the case of the conductor pulls down on the right side. This influential pair is widely separated, affecting the coefficient of income by making it smaller than it would otherwise be. The case of railroad engineers appears as a very high leverage point due to their unusual income given their education level, but it aligns more or less with the trend observed throughout the data. Similarly, the case of the reporter also follows the overall trend.

Similarly, in the education AV plots, we observe a similar trend except for the effect of the minister’s case, which pulls up on the left side, leading to a larger coefficient of education than it would otherwise be. Again, the cases of railroad engineers and reporters seem to align with the general trend of the data, requiring less concern.

Overall, based on the AV plots, we may speculate that the case of the minister could be an influential data point affecting our model’s estimation. We can proceed by updating our original regression model and then comparing the coefficients between the models to validate this hypothesis.

m.duncan2 <- lm(prestige ~ income + education, 
                data = subset(Duncan, !(rownames(Duncan) %in% c("minister"))))

	Model 1			Model 2
Variables	Estimates	S.E.	C.I. (95%)	Estimates	S.E.	C.I. (95%)
(Intercept)	-6.06	4.27	-14.69 – 2.56	-6.63	3.89	-14.48 – 1.22
income	0.60 ***	0.12	0.36 – 0.84	0.73 ***	0.12	0.50 – 0.97
education	0.55 ***	0.10	0.35 – 0.74	0.43 ***	0.10	0.24 – 0.63
Observations	45			44
R2 / R2 adjusted	0.828 / 0.820			0.856 / 0.849
p<0.05   ** p<0.01   *** p<0.001

Indeed, after excluding the case of the minister, the coefficient of income increases from 0.6 to 0.73, while the coefficient of education decreases from 0.55 to 0.43.

We can also run the influencePlot() to check if the Cook’s distance of the minister’s case is indeed very problematic or not.

##                StudRes        Hat      CookD
## minister     3.1345186 0.17305816 0.56637974
## reporter    -2.3970224 0.05439356 0.09898456
## conductor   -1.7040324 0.19454165 0.22364122
## RR.engineer  0.8089221 0.26908963 0.08096807

2. Diagnostics for Linearity and Errors

2.1 Violation of Linearity Assumption

Recall the linearity assumption is that the expected (mean) value of the disturbance term is 0 (i.e., \(E(u_i)=0\)). The average regression error is 0 everywhere implies that the regression surface captures the dependency of the conditional mean of \(y\) on \(x\)s. Violating the linearity assumption implies “that the model fails to capture the systematic pattern of relationship between the response and explanatory variables” (Fox, 2016: 307). For example, a partial relationship specified to be linear may be nonlinear, or two explanatory variables specified to have additive partial effects may interact in determining \(y\). While a fitted model may sometimes provide a useful approximation to the true regression surface \(E(y)\), there are instances where the model can be highly misleading. You can think of nonlinearity (fitting the wrong equation to the data) as potentially the most serious problem with a regression model.

When performing a linear regression with a single explanatory variable, a scatterplot of the response variable against the explanatory variable provides a good indication of the nature of the relationship. If there is more than one explanatory variable, things become more complicated. Although it can still be useful to generate scatterplots of the response variable against each of the explanatory variables, this does not take into account the effect of the other explanatory variables in the model. After all, our interest centers on the partial relationship between \(y\) and each \(x\), controlling for the other \(x\)s, not on the marginal relationship between y and a single \(x\).

Furthermore, while plotting residuals against each \(x\) is helpful for detecting departures from linearity, residual plots cannot distinguish between monotone and nonmonotone nonlinearity.

In the above figure from the Fox textbook, case (a) might be modeled by \(y = \beta_0 + \beta_1\sqrt{x} + \varepsilon\) (a transformation of \(x\)), but case (b) cannot be linearied by a transformation of \(x\) and might instead be dealt with a quadratic regression, \(y = \beta_0 + \beta_1x + \beta_2x^2 + \varepsilon\).

Component + residual (C+R) plots, an extension of partial residual plots, are a good way to see if each predictor you include in the multivariate linear model have a linear relationship to the dependent variable. A partial residual plot essentially attempts to model the residuals of one predictor against the dependent variable. The partial residual for the \(j\)th explanatory variables is defined as:

\[E_{i}^{(j)}= E_i+B_jX_{ij}\]

In words, add back the linear component of the partial relationship between \(Y\) and \(X_j\) to the least-squares residuals, which may include an unmodeled nonlinear component. Then plot \(E^{(j)}\) versus \(X_j\).

Example

Let’s do a simulation example to see how C+R plots can help distinguish monotone and nonmonotone nonlinearity.

set.seed(1234) # for replication
N <- 100 # set up the sample size
x1 <- rnorm(N, mean = 4, sd = 5)
x2 <- runif(N, min = -5, max = 5)
x3 <- rnorm(N, mean = 10, sd = 2)*(x1/2)
cor(x1, x3)

## [1] 0.9724877

e <- rnorm(N, mean = 0, sd = 1) # set the errors to be normally distributed 
y <- 7 + x1^2 + x2 + 5*x2^2 - 7*x3 + e # set the true y
df <- data.frame(y = y, x1 = x1, x2 = x2, x3 = x3, e = e)

In the above chunk, I define the true value \(y\) as a function of \(x_1\), \(x_2\), \(x_3\), and the error term (\(e\)). Specifically, the true regression coefficients are set as follows: \(1\) for the squared transformation of \(x_1\), \(1\) for \(x_2\), and \(-7\) for \(x_3\). Additionally, the term \(5 \cdot x_2^2\) signifies a quadratic effect of \(x_2\), while the intercept is fixed at \(7\). Note, I sepcifically set \(x_3\) to be highly correlated with \(x_1\). To visualize the underlying relationship between \(x_1\), \(x_2\), \(x_3\), and \(y\), let’s start by plotting scatterplots for each pair of \(y\) and \(x\), then running a multivariate regression, and plotting residuals against each \(x\).

fit <- lm(y ~ x1 + x2 + x3, data = df)
df$residual <- fit$residuals

library(tidyverse)

p1 <- ggplot(data = df, aes(x = x1, y = y)) +
  geom_point(color = "gray37", alpha = 0.7) +
  geom_smooth(color = "salmon2", alpha = 0.3) +
  labs(subtitle = "Scatterplot between y and each x") +
  theme_bw()
  
p2 <- ggplot(data = df, aes(x = x1, y = residual)) +
  geom_point(color = "gray37", alpha = 0.7) +
  geom_smooth(method = "lm", color = "steelblue2", alpha = 0.3) +
  labs(subtitle = "Residuals against each x") +
  theme_bw()

p3 <- ggplot(data = df, aes(x = x2, y = y)) +
  geom_point(color = "gray37", alpha = 0.7) +
  geom_smooth(color = "salmon2", alpha = 0.3) +
  theme_bw()
  
p4 <- ggplot(data = df, aes(x = x2, y = residual)) +
  geom_point(color = "gray37", alpha = 0.7) +
  geom_smooth(method = "lm", color = "steelblue2", alpha = 0.3) +
  theme_bw()

p5 <- ggplot(data = df, aes(x = x3, y = y)) +
  geom_point(color = "gray37", alpha = 0.7) +
  geom_smooth(color = "salmon2", alpha = 0.3) +
  theme_bw()
  
p6 <- ggplot(data = df, aes(x = x3, y = residual)) +
  geom_point(color = "gray37", alpha = 0.7) +
  geom_smooth(method = "lm", color = "steelblue2", alpha = 0.3) +
  theme_bw()
  
ggpubr::ggarrange(p1, p2, p3, p4, p5, p6, ncol = 2, nrow = 3)

From the scatterplots of \(y\) against \(x_1\) and \(x_2\), as well as their respective residual plots, it’s evident that while plotting residual plots can indicate whether our model violates the linearity assumption, it doesn’t provide insight into the underlying reasons for the violation. Notably, the residuals plotted against both \(x_1\) and \(x_2\) exhibit a similar convex pattern. Moreover, scatterplots of \(y\) against \(x_1\) and \(x_2\) do not accurately depict their partial relationship. Simply plotting scatterplots between variables fails to consider the influence of other variables in the model, particularly when there are high correlations among predictors. Let’s do the C+R plots instead.

crPlots(fit)

The magenta line represents a smooth curve through the C+R plots, while the blue dashed line represents the regression coefficient of each \(x\) from the multivariate regression model. A notable disparity between the residual line and the component line suggests that the predictor does not exhibit a linear relationship with the response variable. Once again, the advantage of using C+R plots is that they provide information about the partial relationship between \(y\) and each \(x\) while controlling for all other variables. Additionally, they offer insight into how to address nonlinearity. To address the nonlinearity in \(x_1\), applying a square transformation often resolves the issue effectively. As for \(x_2\), we need to consider the quadratic relationship between \(y\) and \(x_2\).

df$x1_square <- (df$x1)^2
fit_update <- lm(y ~ x1_square + poly(x2, 2, raw = TRUE) + x3, data = df)
crPlots(fit_update)

2.2 Violation of Homoskedasticity Assumption

Note: the materials for this section draw from Penn State STAT 501 and here.

The method of ordinary least squares assumes that there is constant variance in the errors (which is called homoscedasticity). The method of weighted least squares can be used when the ordinary least squares assumption of constant variance in the errors is violated (which is called heteroscedasticity). The model under consideration is:

\[Y = X\beta + \epsilon^*\] where \(\epsilon^*\) is assumed to be (multivariate) normally distributed with mean vector 0 and nonconstant variance-covariance matrix.

\[\begin{equation*} \left(\begin{array}{cccc} \sigma^{2}_{1} & 0 & \ldots & 0 \\ 0 & \sigma^{2}_{2} & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \ldots & \sigma^{2}_{n} \\ \end{array} \right) \end{equation*}\]

If we define the reciprocal of each variance, \(\sigma_i^2\), as the weights, \(w_i = 1/\sigma_i^2\), then let matrix W be a diagonal matrix containing these weights:

\[\begin{equation*}\textbf{W}=\left( \begin{array}{cccc} w_{1} & 0 & \ldots & 0 \\ 0& w_{2} & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0& 0 & \ldots & w_{n} \\ \end{array} \right) \end{equation*}\] The weighted least squares estimate is then:

\[\begin{align*} \hat{\beta}_{WLS}&=\arg\min_{\beta}\sum_{i=1}^{n}\epsilon_{i}^{*2}\\ &=(\textbf{X}^{T}\textbf{W}\textbf{X})^{-1}\textbf{X}^{T}\textbf{W}\textbf{Y} \end{align*}\] With this setting, we can make a few observations:

• Since each weight is inversely proportional to the error variance, it reflects the information in that observation. So, an observation with a small error variance has a large weight since it contains relatively more information than an observation with a large error variance (small weight).
• The weights have to be known (or more usually estimated) up to a proportionality constant.
• Weighted least squares estimates of the coefficients will usually be nearly the same as the “ordinary” unweighted estimates. In cases where they differ substantially, the procedure can be iterated until estimated coefficients stabilize (often in no more than one or two iterations); this is called iteratively reweighted least squares.
• The use of weights will (legitimately) impact the widths of statistical intervals.

Example

Let’s simulate the violation of homoskedasticity assumption and apply the method of weighted least squares to address it.

set.seed(1234) # for replication
N <- 1000 # set up the sample size
x1 <- 1:N
x2 <- rnorm(N, mean = 10, sd = 3)
sd <- runif(N, min = 0, max = 2)
e <- rnorm(N, mean = 0, sd = 2*x1) # set the errors to be normally distributed but not at a constant variance
y <- 7 + 3*x1 + x2 + e # set the true y
df <- data.frame(y = y, x1 = x1, x2 = x2, e = e)

In the above chunk, I set the error variance to be more pronounce: as \(x1\) gets larger, the error variance gets larger by multiplying \(x1\). You can see this relationship in the below figure.

ggplot(data = df, aes(x = x1, y = e)) +
  geom_point(color = "gray37", alpha = 0.7) +
  theme_bw()

Let’s now run the regression and test it to see if it violates the homoskedasticity assumption.

fit <- lm(y ~ x1 + x2, data = df) 

	OLS
Variables	Estimates	S.E.	C.I. (95%)
(Intercept)	-65.82	144.57	-349.52 – 217.87
x1	3.13 ***	0.13	2.88 – 3.38
x2	6.73	12.46	-17.71 – 31.17
Observations	1000
R2 / R2 adjusted	0.371 / 0.370
p<0.05   ** p<0.01   *** p<0.001

According to the regression table, you can see that the standard error of \(x2\) is very large, which is due to the heteroskedasticity problem in our data. The estimate of the effect of \(x2\), therefore, is very inefficient.

par(mfrow = c(2, 2))
plot(fit)

Let’s take a look at the third plot. Scale-Location is used to check the homogeneity of variance of the residuals (homoscedasticity). Horizontal line with equally spread points is a good indication of homoscedasticity. This is not the case in our example, where we have a heteroscedasticity problem. It can be seen that the variability (variances) of the residual points increases with the value of the fitted outcome variable, suggesting non-constant variances in the residuals errors.

To formally test for heteroscedasticity, we can perform a Breusch-Pagan test:

# load lmtest package
library(lmtest)

# perform Breusch-Pagan test
bptest(fit)

## 
##  studentized Breusch-Pagan test
## 
## data:  fit
## BP = 168.84, df = 2, p-value < 2.2e-16

The Breusch-Pagan test uses the following null and alternative hypotheses:

• Null Hypothesis (\(H_0\)): Homoscedasticity is present (the residuals are distributed with equal variance)
• Alternative Hypothesis (\(H_A\)): Heteroscedasticity is present (the residuals are not distributed with equal variance)

Since the p-value from the test is < 2.2e-16, we will reject the null hypothesis and conclude that heteroscedasticity is a problem in this model.

Since heteroscedasticity is present, we will perform weighted least squares by defining the weights in such a way that the observations with lower variance are given more weight:

# define weights to use
wt <- 1 / lm(abs(fit$residuals) ~ fit$fitted.values)$fitted.values^2

head(lm(abs(fit$residuals) ~ fit$fitted.values)$fitted.values)

##        1        2        3        4        5        6 
## 22.48226 39.03961 48.75519 15.68057 45.26871 47.61134

head(wt)

##            1            2            3            4            5            6 
## 0.0019784279 0.0006561287 0.0004206862 0.0040670185 0.0004879820 0.0004411429

# perform weighted least squares regression
wls_fit <- lm(y ~ x1 + x2, data = df, weights = wt)

	No WLS			WLS
Variables	Estimates	S.E.	C.I. (95%)	Estimates	S.E.	C.I. (95%)
(Intercept)	-65.82	144.57	-349.52 – 217.87	2.14	20.99	-39.04 – 43.32
x1	3.13 ***	0.13	2.88 – 3.38	3.09 ***	0.07	2.94 – 3.23
x2	6.73	12.46	-17.71 – 31.17	0.79	2.55	-4.21 – 5.79
Observations	1000			1000
R2 / R2 adjusted	0.371 / 0.370			0.648 / 0.648
p<0.05   ** p<0.01   *** p<0.001

The weighted least squares model yields a coefficient for \(x_2\) much closer to its true value compared to the original model. Additionally, the residual standard error in the WLS model is smaller than that in the original multivariate linear regression model. These results suggest that the predicted values generated by the weighted least squares model are much closer to the actual observations than those produced by the original model.

The weighted least squares model also has an R-squared of 0.648 compared to 0.371 in the original multivariate linear regression model. This indicates that the weighted least squares model is able to explain more of the variance in y compared to the original multivariate linear regression model.

Discussion IV

For today’s discussion, we will review what multicollinearity is, how to detect it, and introduce Principal Component Analysis (PCA) as one of the methods to handle it.

1. (Multi)collinearity Problem
2. Detecting Collinearity
3. Principle Component Analysis (PCA)

1. (Multi)collinearity Problem

In multivariate linear regression, it is possible for two or more predictor variables to be correlated with each other. This scenario is referred to as collinearity. In extreme cases, where collinearity exists between three or more variables, even if no pair of variables has a particularly high correlation, this is known as multicollinearity. Multicollinearity indicates redundancy between predictor variables.

The presence of multicollinearity can cause the solution of the regression model to become unstable.

• coefficient standard errors are large reflecting imprecision
• broad confidence intervals
• hypothesis test have low power

Therefore, it is important to detect and address multicollinearity before building a regression model to ensure the model’s reliability and accuracy.

1.1 Example of Perfect Multicollinearity Problem

Let’s first review what Lauren discussed in the lecture regarding why it’s not feasible to include all the dummies of a categorical variable .

library(AER)
library(MASS)
data(CASchools)

# define variables
CASchools$STR <- CASchools$students/CASchools$teachers 
CASchools$score <- (CASchools$read + CASchools$math)/2 

# set seed for reproducibility
set.seed(1)

# generate artificial data on location
CASchools$direction <- sample(c("West", "North", "South", "East"), 420, replace = T)
head(CASchools[, c("score", "STR", "english", "direction")])

##    score      STR   english direction
## 1 690.80 17.88991  0.000000      West
## 2 661.20 21.52466  4.583333      East
## 3 643.60 18.69723 30.000002     South
## 4 647.70 17.35714  0.000000      West
## 5 640.85 18.67133 13.857677     North
## 6 605.55 21.40625 12.408759      West

In this example, our response variable is score, and explanatory variables are STR(student-teacher ratio), english (the share of English learning students), and direction (artificial location of a school).

mod1 <- lm(score ~ STR + english + direction, data = CASchools)

	Score
Variables	Estimates	S.E.	C.I. (95%)
(Intercept)	684.80 ***	7.54	669.98 – 699.63
STR	-1.09 **	0.38	-1.84 – -0.34
english	-0.66 ***	0.04	-0.73 – -0.58
direction [North]	1.66	2.06	-2.38 – 5.71
direction [South]	0.72	2.06	-3.34 – 4.77
direction [West]	1.79	1.98	-2.10 – 5.69
Observations	420
R2 / R2 adjusted	0.428 / 0.421
p<0.05   ** p<0.01   *** p<0.001

In the regression table above, you’ll notice that R has included only three of the dummies for the direction variable and has designated the category East as the reference group in the model. What this means is that when interpreting the coefficients of the remaining three categories of the direction variable, we are comparing their effects to the category East. For instance, we can infer that schools located in the North have, on average, 1.66 higher test scores than schools located in the East. Although this difference is not statistically significant, as the variable was artificially generated.

In addition, since the category East was designated as the reference group, the value of the intercept reflects that the baseline test score of a school is when schools have a student-teacher ratio of 0, no English learning students, and are located in the East. The modelpredicts that the average test score of such schools is 684.80.

What would happen if we intentionally include four dummies of the direction variable in the model?

library(fastDummies)
CASchools <- dummy_cols(CASchools, select_columns = "direction")
head(CASchools[,c("score", "STR", "english", "direction_East", 
                  "direction_North", "direction_South", "direction_West")])

##    score      STR   english direction_East direction_North direction_South
## 1 690.80 17.88991  0.000000              0               0               0
## 2 661.20 21.52466  4.583333              1               0               0
## 3 643.60 18.69723 30.000002              0               0               1
## 4 647.70 17.35714  0.000000              0               0               0
## 5 640.85 18.67133 13.857677              0               1               0
## 6 605.55 21.40625 12.408759              0               0               0
##   direction_West
## 1              1
## 2              0
## 3              0
## 4              1
## 5              0
## 6              1

mod2 <- lm(score ~ STR + english + direction_North + direction_South + 
               direction_West + direction_East , data = CASchools)

	Model 1			Model 2
Variables	Estimates	S.E.	C.I. (95%)	Estimates	S.E.	C.I. (95%)
(Intercept)	684.80 ***	7.54	669.98 – 699.63	684.80 ***	7.54	669.98 – 699.63
STR	-1.09 **	0.38	-1.84 – -0.34	-1.09 **	0.38	-1.84 – -0.34
english	-0.66 ***	0.04	-0.73 – -0.58	-0.66 ***	0.04	-0.73 – -0.58
direction [North]	1.66	2.06	-2.38 – 5.71
direction [South]	0.72	2.06	-3.34 – 4.77
direction [West]	1.79	1.98	-2.10 – 5.69
direction North				1.66	2.06	-2.38 – 5.71
direction South				0.72	2.06	-3.34 – 4.77
direction West				1.79	1.98	-2.10 – 5.69
Observations	420			420
R2 / R2 adjusted	0.428 / 0.421			0.428 / 0.421
p<0.05   ** p<0.01   *** p<0.001

As you can observe in model 2, R automatically dropped one of the dummies of the direction variable (East). R does this to solve the perfect multicollinearity problem. Where is the perfect multicollinearity problem in this case?

Recall from the materials in Discussion II, we talk about how to think of our data set as a matrix. The below equations are how each observation in our data set generated.

\(Y_1=\beta_0+\beta_1X_{11}+\beta_2X_{12}+...+\beta_kX_{1k}+\varepsilon_1\)

\(Y_2=\beta_0+\beta_1X_{21}+\beta_2X_{22}+...+\beta_kX_{2k}+\varepsilon_2\)

.

.

.

\(Y_n=\beta_0+\beta_1X_{n1}+\beta_2X_{n2}+...+\beta_kX_{nk}+\varepsilon_n\)

Rewrite them into a matrix form.

The system of equations is summarized as the y vector equal to the product of the X matrix of data times the \(\beta\) vector of parameters plus the vector of disturbances, \(\varepsilon\). In the lecture, Lauren showed how to apply what we did when finding a pair of intercept and slope in simple linear regression to solve the optimization problem to find the coefficient vector b using matrix algebra. Instead of taking derivative with respect to each coefficient one at a time, we can take derivative directly with respect to the coefficient vector b. Therefore, the least squares estimator is:

\(\hat{\beta}=(X'X)^{-1}X'y\)

In the matrix algebra for estimating \(\beta\) coefficients using the above equation, it is necessary for the first column of the \(X\) matrix to be composed of 1s. This is because the first column represents the intercept term in the regression model. Therefore, including all the dummies of a categorical variable in a regression model would result in a perfect multicollinearity problem beacuse they are the linear combination of the intercept term.

Let’s relook at the CASchools again.

head(CASchools[,c("score", "STR", "english", "direction_East", 
                  "direction_North", "direction_South", "direction_West")])

##    score      STR   english direction_East direction_North direction_South
## 1 690.80 17.88991  0.000000              0               0               0
## 2 661.20 21.52466  4.583333              1               0               0
## 3 643.60 18.69723 30.000002              0               0               1
## 4 647.70 17.35714  0.000000              0               0               0
## 5 640.85 18.67133 13.857677              0               1               0
## 6 605.55 21.40625 12.408759              0               0               0
##   direction_West
## 1              1
## 2              0
## 3              0
## 4              1
## 5              0
## 6              1

The linear combination of the intercept term equals to the sum of direction_East, direction_North, direction_South, direction_West, resulting in a perfect multicollinearity problem.

CASchools$intercept <- CASchools$direction_East + CASchools$direction_North + 
    CASchools$direction_South + CASchools$direction_West

head(CASchools[,c("score", "STR", "english", "direction_East", 
                  "direction_North", "direction_South", "direction_West", "intercept")])

##    score      STR   english direction_East direction_North direction_South
## 1 690.80 17.88991  0.000000              0               0               0
## 2 661.20 21.52466  4.583333              1               0               0
## 3 643.60 18.69723 30.000002              0               0               1
## 4 647.70 17.35714  0.000000              0               0               0
## 5 640.85 18.67133 13.857677              0               1               0
## 6 605.55 21.40625 12.408759              0               0               0
##   direction_West intercept
## 1              1         1
## 2              0         1
## 3              0         1
## 4              1         1
## 5              0         1
## 6              1         1

1.2 Example of Imperfect Multicollinearity Problem

If there is high multicollinearity between variables but not perfectly correlated, lm() will still spit out the estimates of the highly correlated variable. However, the standard errors of coefficients will be very high, which means that our estimates will be highly unstable.

set.seed(1234)

# generate a highly correlated variable
CASchools$english_hc1 <- 2*CASchools$english + rnorm(nrow(CASchools), 10, 5)

cor(CASchools$english_hc1, CASchools$english) # english_hc is highly but not perfectly correlated with english

## [1] 0.9904646

set.seed(1234)

# generate a less highly correlated variable
CASchools$english_hc2 <- 2*CASchools$english + rnorm(nrow(CASchools), 10, 20)

cor(CASchools$english_hc2, CASchools$english) 

## [1] 0.8662955

mod3 <- lm(score ~ STR + english , data = CASchools) # without highly correlated term
mod4 <- lm(score ~ STR + english + english_hc1 , data = CASchools) # with highly correlated term
mod5 <- lm(score ~ STR + english + english_hc2 , data = CASchools) # with less highly correlated term

	Model 1			Model 2			Model 3
Variables	Estimates	S.E.	C.I. (95%)	Estimates	S.E.	C.I. (95%)	Estimates	S.E.	C.I. (95%)
(Intercept)	686.03 ***	7.41	671.46 – 700.60	687.24 ***	7.57	672.37 – 702.12	686.40 ***	7.43	671.80 – 701.00
STR	-1.10 **	0.38	-1.85 – -0.35	-1.10 **	0.38	-1.85 – -0.36	-1.10 **	0.38	-1.85 – -0.36
english	-0.65 ***	0.04	-0.73 – -0.57	-0.43	0.28	-0.98 – 0.12	-0.60 ***	0.08	-0.75 – -0.44
english hc1				-0.11	0.14	-0.39 – 0.16
english hc2							-0.03	0.04	-0.10 – 0.04
Observations	420			420			420
R2 / R2 adjusted	0.426 / 0.424			0.427 / 0.423			0.427 / 0.423
p<0.05   ** p<0.01   *** p<0.001

The regression output indicates that in models 2 and 3, there is an imperfect multicollinearity problem, which is reflected by an increase in the standard error of the coefficient estimates. Specifically, the estimated effect of english becomes increasingly unstable as the correlation between english and the correlated term (which we artificially generated) becomes stronger.

2. Detecting Collinearity

Multicollinearity in a regression model can be assessed by computing the variance inflation factor (VIF) for each predictor variable. VIF measures how much of the variation in one variable is explained by the other variable. This is done by running a regression using one of the correlated \(x\) variables as the dependent variable against the other variables as predictor variables. The VIF is calculated as one divided by the tolerance, which is defined as one minus R-squared.

\(VIF_j = \frac{1}{1-R_j^2}\)

A VIF value of 1 indicates an absence of multicollinearity, while values exceeding 5 or 10 suggest problematic collinearity (James et al., 2014). When multicollinearity is detected, it is recommended to remove one or more of the correlated variables from the model. This is because the presence of multicollinearity implies that the information provided by the correlated variables about the response variable is redundant in the presence of other variables (James et al., 2014; P. Bruce & Bruce, 2017).

The R function vif() in the car package can be used to detect multicollinearity in a regression model:

library(car)
vif(mod3) # without highly correlated term

##      STR  english 
## 1.036495 1.036495

vif(mod4) # with highly correlated term

##         STR     english english_hc1 
##    1.036565   52.750311   52.691167

vif(mod5) # with less highly correlated term

##         STR     english english_hc2 
##    1.036565    4.049684    4.007776

As anticipated, when we include a highly correlated term in the model, both english and english_hc1 exhibit extremely high VIF scores.

3. Principal Components Analysis (PCA)

Principal Components Analysis (PCA), which is an unsupervised machine learning technique, can help handle multicollinearity problems by creating new features that are orthogonal to each other, meaning that they have zero correlation. This eliminates the problem of having features that are redundant or that inflate the variance of your model coefficients.

The goal of PCA is to explain most of the variability in a dataset with fewer variables than the original dataset. PCA offers a way to find a low-dimensional representation of a dataset that captures as much of the variation in the data as possible. If we’re able to capture most of the variation in just two dimensions, we could project all of the observations in the original dataset onto a simple scatterplot.

The way we find the principal components is as follows:

Given a dataset with \(p\) predictors: \(X_1\), \(X_2\),…, \(X_p\), calculate \(Z_1\),…, \(Z_M\) to be the \(M\) linear combinations of the original \(p\) predictors where:

• \(Z_m = \Sigma \Phi_{jm}X_j\) for some constants \(\Phi_{1m}\), \(\Phi_{2m}\),…, \(\Phi_{pm}\), \(m = 1, 2,...,M\).
• \(Z_1\) is the linear combination of the predictors that captures the most variance possible.
• \(Z_2\) is the next linear combination of the predictors that captures the most variance while being orthogonal (i.e. uncorrelated) to \(Z_1\).
• \(Z_3\) is then the next linear combination of the predictors that captures the most variance while being orthogonal to \(Z_2\).
• And so on.

In practice, we use the following steps to calculate the linear combinations of the original predictors:

1. Scale each of the variables to have a mean of 0 and a standard deviation of 1.

2. Calculate the covariance matrix for the scaled variables.

3. Calculate the eigenvalues of the covariance matrix.

Using linear algebra, it can be shown that the eigenvector that corresponds to the largest eigenvalue is the first principal component. In other words, this particular combination of the predictors explains the most variance in the data.

The eigenvector corresponding to the second largest eigenvalue is the second principal component, and so on.

Example

For this example, we will use the “Communities and Crime” dataset from the UCI Machine Learning Repository, created by Michael Redmond. This dataset includes various socio-economic, law enforcement, and crime data. A cleaned version is available in the mogavs package under the name crimeData, which has been simplified by removing all columns identified as redundant by the dataset’s author, reducing the number of columns from 128 to 123.

We will use this dataset from the mogavs package for our analysis. The outcome variable we are focusing on is violent crimes per capita, while the independent variables include aspects of law enforcement, socio-economic factors, and crime data. A detailed description of all these variables can be found here.

library(mogavs)
dim(crimeData)

## [1] 1994  123

colnames(crimeData)

##   [1] "x.V6"   "x.V7"   "x.V8"   "x.V9"   "x.V10"  "x.V11"  "x.V12"  "x.V13" 
##   [9] "x.V14"  "x.V15"  "x.V16"  "x.V17"  "x.V18"  "x.V19"  "x.V20"  "x.V21" 
##  [17] "x.V22"  "x.V23"  "x.V24"  "x.V25"  "x.V26"  "x.V27"  "x.V28"  "x.V29" 
##  [25] "x.V30"  "x.V31"  "x.V32"  "x.V33"  "x.V34"  "x.V35"  "x.V36"  "x.V37" 
##  [33] "x.V38"  "x.V39"  "x.V40"  "x.V41"  "x.V42"  "x.V43"  "x.V44"  "x.V45" 
##  [41] "x.V46"  "x.V47"  "x.V48"  "x.V49"  "x.V50"  "x.V51"  "x.V52"  "x.V53" 
##  [49] "x.V54"  "x.V55"  "x.V56"  "x.V57"  "x.V58"  "x.V59"  "x.V60"  "x.V61" 
##  [57] "x.V62"  "x.V63"  "x.V64"  "x.V65"  "x.V66"  "x.V67"  "x.V68"  "x.V69" 
##  [65] "x.V70"  "x.V71"  "x.V72"  "x.V73"  "x.V74"  "x.V75"  "x.V76"  "x.V77" 
##  [73] "x.V78"  "x.V79"  "x.V80"  "x.V81"  "x.V82"  "x.V83"  "x.V84"  "x.V85" 
##  [81] "x.V86"  "x.V87"  "x.V88"  "x.V89"  "x.V90"  "x.V91"  "x.V92"  "x.V93" 
##  [89] "x.V94"  "x.V95"  "x.V96"  "x.V97"  "x.V98"  "x.V99"  "x.V100" "x.V101"
##  [97] "x.V102" "x.V103" "x.V104" "x.V105" "x.V106" "x.V107" "x.V108" "x.V109"
## [105] "x.V110" "x.V111" "x.V112" "x.V113" "x.V114" "x.V115" "x.V116" "x.V117"
## [113] "x.V118" "x.V119" "x.V120" "x.V121" "x.V122" "x.V123" "x.V124" "x.V125"
## [121] "x.V126" "x.V127" "y"

# initialize
crimeData.y <- crimeData$y # set y as a vector out of the original data  
crimeData.X <- crimeData[, -123] # takeout y

By plotting the correlation matrix among predictors, we observe significant correlations between some variables, suggesting potential multicollinearity issues if all variables are included in the model.

res <- cor(crimeData.X, method = "pearson")
corrplot::corrplot(res, method= "color", 
                   order = "hclust", tl.pos = 'n')

As you compute the Variance Inflation Factor (VIF) for each predictor, you will notice that many variables have VIF scores exceeding 5 or 10, indicating potentially problematic multicollinearity.

fit <- lm(y ~ ., data = crimeData)
summary(fit)$r.squared

## [1] 0.6979036

round(vif(fit), 2)

##      x.V6      x.V7      x.V8      x.V9     x.V10     x.V11     x.V12     x.V13 
##    411.21     50.92     27.09     43.71      8.04     27.76     56.63     80.66 
##     x.V14     x.V15     x.V16     x.V17     x.V18     x.V19     x.V20     x.V21 
##    164.73     76.54    413.92     20.03    237.53     66.08      4.75     24.54 
##     x.V22     x.V23     x.V24     x.V25     x.V26     x.V27     x.V28     x.V29 
##     61.59     18.88      7.07    138.42    187.64    107.17      7.55      1.83 
##     x.V30     x.V31     x.V32     x.V33     x.V34     x.V35     x.V36     x.V37 
##      2.69      2.94      2.79     40.68     42.25     32.56     59.08     52.65 
##     x.V38     x.V39     x.V40     x.V41     x.V42     x.V43     x.V44     x.V45 
##     13.33     41.26      5.93      8.64     19.69     44.51    295.17     27.51 
##     x.V46     x.V47     x.V48     x.V49     x.V50     x.V51     x.V52     x.V53 
##    442.21   1319.59    111.72    166.85    165.81     22.10     16.23     11.73 
##     x.V54     x.V55     x.V56     x.V57     x.V58     x.V59     x.V60     x.V61 
##     14.03     18.61     16.77      6.46     23.16     54.16     64.93     32.00 
##     x.V62     x.V63     x.V64     x.V65     x.V66     x.V67     x.V68     x.V69 
##    137.85    460.77    732.28    465.34     45.15     37.96    285.81    302.77 
##     x.V70     x.V71     x.V72     x.V73     x.V74     x.V75     x.V76     x.V77 
##    306.56    116.33     49.29    732.40     40.51     17.87      3.60     14.91 
##     x.V78     x.V79     x.V80     x.V81     x.V82     x.V83     x.V84     x.V85 
##      5.69    711.20      3.02      3.90      6.24     11.52      3.67    347.53 
##     x.V86     x.V87     x.V88     x.V89     x.V90     x.V91     x.V92     x.V93 
##    777.26    203.14     37.64    249.43     82.63    133.85      5.08      6.59 
##     x.V94     x.V95     x.V96     x.V97     x.V98     x.V99    x.V100    x.V101 
##      3.50      5.35      2.71     75.55     12.98     19.82     11.48     12.03 
##    x.V102    x.V103    x.V104    x.V105    x.V106    x.V107    x.V108    x.V109 
##    292.31 558345.67    175.45    438.06     16.57     45.50     14.09 552634.18 
##    x.V110    x.V111    x.V112    x.V113    x.V114    x.V115    x.V116    x.V117 
##     12.67     14.35     55.13     51.63      8.50    106.64     16.77      3.99 
##    x.V118    x.V119    x.V120    x.V121    x.V122    x.V123    x.V124    x.V125 
##      4.57      3.63      5.15      3.95      9.26     51.25     19.37      3.65 
##    x.V126    x.V127 
##      2.67     68.40

Instead of arbitrarily dropping variables from our model to address multicollinearity, we can perform principal component analysis (PCA). This approach allows us to identify principal components that capture most of the variance in the data with fewer variables, thus avoiding multicollinearity and leading to more efficient estimates.

Now, the first procedure for PCA is to scale and center the data using prcomp command. This command is both convenient and useful, helping you to subtract each observation value by its mean and then divide the result by the standard deviation. Then, we can use print and summary commands to see the rotation, standard deviation and proportion of variance each component contributes to the dataset, and use fviz_fig in package factoextra to plot a nice formatted scree plot.

pca <- prcomp(crimeData.X, center = TRUE, scale. = TRUE)
summary(pca)

## Importance of components:
##                           PC1    PC2     PC3    PC4     PC5     PC6    PC7
## Standard deviation     5.1263 4.2852 3.14023 2.8484 2.52280 2.44983 2.1274
## Proportion of Variance 0.2154 0.1505 0.08083 0.0665 0.05217 0.04919 0.0371
## Cumulative Proportion  0.2154 0.3659 0.44674 0.5132 0.56542 0.61461 0.6517
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     1.92926 1.79021 1.56346 1.46739 1.39980 1.30804 1.28740
## Proportion of Variance 0.03051 0.02627 0.02004 0.01765 0.01606 0.01402 0.01359
## Cumulative Proportion  0.68222 0.70849 0.72852 0.74617 0.76223 0.77626 0.78984
##                           PC15    PC16    PC17    PC18    PC19    PC20    PC21
## Standard deviation     1.26287 1.19862 1.14023 1.12564 1.05175 0.99584 0.97054
## Proportion of Variance 0.01307 0.01178 0.01066 0.01039 0.00907 0.00813 0.00772
## Cumulative Proportion  0.80291 0.81469 0.82535 0.83573 0.84480 0.85293 0.86065
##                           PC22    PC23    PC24    PC25    PC26    PC27    PC28
## Standard deviation     0.96633 0.92222 0.89262 0.85634 0.84920 0.82465 0.77804
## Proportion of Variance 0.00765 0.00697 0.00653 0.00601 0.00591 0.00557 0.00496
## Cumulative Proportion  0.86830 0.87527 0.88181 0.88782 0.89373 0.89930 0.90426
##                           PC29    PC30    PC31    PC32   PC33    PC34    PC35
## Standard deviation     0.76580 0.73427 0.72819 0.70541 0.6899 0.68415 0.67038
## Proportion of Variance 0.00481 0.00442 0.00435 0.00408 0.0039 0.00384 0.00368
## Cumulative Proportion  0.90907 0.91349 0.91784 0.92191 0.9258 0.92965 0.93334
##                           PC36    PC37    PC38    PC39    PC40    PC41    PC42
## Standard deviation     0.66397 0.62569 0.60702 0.58972 0.57603 0.56650 0.55393
## Proportion of Variance 0.00361 0.00321 0.00302 0.00285 0.00272 0.00263 0.00252
## Cumulative Proportion  0.93695 0.94016 0.94318 0.94603 0.94875 0.95138 0.95389
##                           PC43    PC44   PC45    PC46    PC47   PC48    PC49
## Standard deviation     0.54304 0.52869 0.5180 0.49787 0.49653 0.4687 0.45403
## Proportion of Variance 0.00242 0.00229 0.0022 0.00203 0.00202 0.0018 0.00169
## Cumulative Proportion  0.95631 0.95860 0.9608 0.96283 0.96486 0.9667 0.96835
##                           PC50    PC51    PC52    PC53   PC54    PC55   PC56
## Standard deviation     0.45164 0.44451 0.43885 0.42861 0.4130 0.41005 0.3987
## Proportion of Variance 0.00167 0.00162 0.00158 0.00151 0.0014 0.00138 0.0013
## Cumulative Proportion  0.97002 0.97164 0.97322 0.97472 0.9761 0.97750 0.9788
##                           PC57    PC58    PC59    PC60   PC61    PC62    PC63
## Standard deviation     0.38490 0.37296 0.37131 0.36040 0.3487 0.32963 0.32066
## Proportion of Variance 0.00121 0.00114 0.00113 0.00106 0.0010 0.00089 0.00084
## Cumulative Proportion  0.98001 0.98115 0.98228 0.98335 0.9843 0.98524 0.98608
##                           PC64    PC65    PC66   PC67    PC68    PC69    PC70
## Standard deviation     0.31633 0.30719 0.29810 0.2913 0.28788 0.27641 0.26830
## Proportion of Variance 0.00082 0.00077 0.00073 0.0007 0.00068 0.00063 0.00059
## Cumulative Proportion  0.98690 0.98767 0.98840 0.9891 0.98978 0.99040 0.99099
##                           PC71    PC72    PC73    PC74    PC75    PC76    PC77
## Standard deviation     0.25982 0.25450 0.25103 0.23801 0.23588 0.22980 0.22319
## Proportion of Variance 0.00055 0.00053 0.00052 0.00046 0.00046 0.00043 0.00041
## Cumulative Proportion  0.99155 0.99208 0.99259 0.99306 0.99351 0.99395 0.99435
##                           PC78    PC79    PC80    PC81    PC82    PC83    PC84
## Standard deviation     0.21832 0.21411 0.21056 0.20020 0.19450 0.18859 0.18256
## Proportion of Variance 0.00039 0.00038 0.00036 0.00033 0.00031 0.00029 0.00027
## Cumulative Proportion  0.99475 0.99512 0.99548 0.99581 0.99612 0.99641 0.99669
##                           PC85    PC86    PC87    PC88   PC89   PC90    PC91
## Standard deviation     0.17283 0.16885 0.16837 0.16166 0.1563 0.1556 0.14596
## Proportion of Variance 0.00024 0.00023 0.00023 0.00021 0.0002 0.0002 0.00017
## Cumulative Proportion  0.99693 0.99717 0.99740 0.99761 0.9978 0.9980 0.99819
##                           PC92    PC93    PC94    PC95    PC96    PC97    PC98
## Standard deviation     0.14342 0.13589 0.13405 0.13088 0.12553 0.11681 0.11606
## Proportion of Variance 0.00017 0.00015 0.00015 0.00014 0.00013 0.00011 0.00011
## Cumulative Proportion  0.99835 0.99851 0.99865 0.99879 0.99892 0.99903 0.99915
##                           PC99   PC100   PC101   PC102   PC103   PC104   PC105
## Standard deviation     0.11383 0.10552 0.10317 0.09484 0.09310 0.08754 0.08136
## Proportion of Variance 0.00011 0.00009 0.00009 0.00007 0.00007 0.00006 0.00005
## Cumulative Proportion  0.99925 0.99934 0.99943 0.99950 0.99957 0.99964 0.99969
##                          PC106   PC107   PC108   PC109   PC110   PC111   PC112
## Standard deviation     0.07377 0.07086 0.06453 0.06077 0.05462 0.05390 0.04860
## Proportion of Variance 0.00004 0.00004 0.00003 0.00003 0.00002 0.00002 0.00002
## Cumulative Proportion  0.99974 0.99978 0.99981 0.99984 0.99987 0.99989 0.99991
##                          PC113   PC114   PC115   PC116   PC117   PC118   PC119
## Standard deviation     0.04829 0.04163 0.03993 0.03839 0.03532 0.02948 0.02670
## Proportion of Variance 0.00002 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001
## Cumulative Proportion  0.99993 0.99994 0.99996 0.99997 0.99998 0.99999 0.99999
##                          PC120   PC121     PC122
## Standard deviation     0.02461 0.02171 0.0009488
## Proportion of Variance 0.00000 0.00000 0.0000000
## Cumulative Proportion  1.00000 1.00000 1.0000000

In the summary part, the line Proportion of Variance represents how much percentage of the variance is explained by each component. Here the first component explains 21.54% of the variability, the second component explains 15.05% of the variance and so on… Then the line Cumulative Proportion represents the percentage of variance explained by first k elements. Specifically, PC41 column gives the value 0.95138, which means if we include the first 41 components, then 95.14% of the variance is explained.

Note that the principal components scores for each observation are stored in pca$x.

# display principal components
head(pca$x)

##            PC1        PC2         PC3        PC4        PC5       PC6
## [1,]  1.522743 -0.7916743 -3.28026359 -1.2753142  0.5157339  2.246354
## [2,] -1.015422  1.0858176 -3.62568762 -3.4731108  3.0243619 -0.354592
## [3,] -3.158742 -2.4840295  0.07796813  1.0533924 -1.5569032 -1.452501
## [4,]  2.289034  2.8487601  2.05171050  3.4190377 -2.7585249  3.700239
## [5,]  5.542577 -2.4608213  1.99899950  1.9062590 -1.0345593  2.313514
## [6,]  6.281678  7.1365764 -3.36227772 -0.4295065  2.8784476 -3.242882
##             PC7        PC8         PC9      PC10      PC11        PC12
## [1,] -2.7985542  0.9267431  0.58399391 0.6004072 0.9551796 -0.79934712
## [2,] -4.4201234  1.9166227 -1.72469679 0.1042599 0.1723816  0.04400804
## [3,] -0.3855891  3.8352238 -0.29908075 2.4358470 2.3923717  1.29196186
## [4,]  2.8111671  5.9181406 -0.77918858 2.7291999 0.3761453  4.26473995
## [5,] -0.6231897 -0.2395864  0.39365486 0.2781169 0.3233575  1.17753888
## [6,]  0.1920078  0.4265712  0.09238477 0.1389794 1.7694596  0.30840918
##            PC13       PC14       PC15       PC16        PC17        PC18
## [1,] -0.9666304  0.1303065 -0.4657864  1.2189250  0.06043401 0.729760127
## [2,] -2.8716352 -0.5208606  1.0407432  0.5515834  1.16689692 0.002392092
## [3,]  0.9311878 -0.8210312  2.1994493  0.7673565 -1.14019003 0.424345038
## [4,]  1.6327709 -2.8870481 -3.5170495  1.0626643 -1.28514836 1.244770149
## [5,] -0.1779068 -0.9654716 -0.2833528 -0.7469055 -0.47950830 0.936034210
## [6,] -2.9758081 -0.3608440  2.0664083  0.5216004  1.48223391 2.949611589
##             PC19       PC20       PC21        PC22          PC23       PC24
## [1,] -0.79921663  0.1319244  0.1489032  0.07092186 -0.3082014815 -0.8177185
## [2,] -0.62406584  0.1334736 -0.3755883  0.12341205  0.9769526455  0.6355450
## [3,] -0.07681935  0.7818357  0.0991699 -0.77126234  0.0405022795  0.2899994
## [4,] -1.60374672  0.9227981  0.1087099  0.35152183  0.0005690714  0.9002704
## [5,]  0.88507386  1.8341636 -0.9523725  0.57316432  1.3744618673  0.8042577
## [6,] -0.70898935 -0.6434944  0.1147150 -0.66398971  0.0099635183  0.5997810
##            PC25       PC26       PC27        PC28        PC29       PC30
## [1,] -0.3069298 -0.4152735  0.2167211  0.21327767 -0.14586468 -0.3158446
## [2,]  0.2369439  1.0503642  0.6764152  0.08486815  0.94130197  1.5350282
## [3,]  0.6215431 -1.3418481  1.8594000  0.31524477 -0.44893741 -0.1837356
## [4,]  0.4220371  1.2271946  1.2904089 -0.67079238  0.08942445 -0.5654088
## [5,] -0.4217490  1.1109450 -1.1807146  0.38248434 -0.11894801  0.6922435
## [6,]  1.2264617 -0.2096174 -0.5534602 -0.04476907  1.08860753  0.2460175
##            PC31        PC32       PC33       PC34        PC35       PC36
## [1,]  0.6630573 -1.02344207  0.3339815 -0.0757706  0.06264455  0.3350997
## [2,]  0.4915675  0.11530529 -0.2695809  0.2543846  0.38000284  0.6615344
## [3,]  0.8941770 -0.20630166 -0.2529353  0.2248543  0.68720542  0.2818908
## [4,]  1.0092294  0.87614189 -0.5055810 -0.6989796  1.60577512 -0.1289466
## [5,]  0.3785551  0.81712031 -0.2265443 -0.4916125  0.59171867 -0.2505556
## [6,] -0.0408304  0.09808626  1.0147548  0.3374860 -0.37652980 -0.4025329
##              PC37       PC38          PC39       PC40        PC41        PC42
## [1,]  0.259200542  0.2418568 -0.0004330426 -0.4746566 -0.14509983 -0.08296911
## [2,]  1.878980303  0.2291362 -0.2598100237 -0.1539732  0.74360461 -0.64468912
## [3,] -0.004971392  0.3541479 -0.4798498803  0.5575868  0.37614191  0.32161952
## [4,]  0.616734544  0.1541563  0.6113958888  0.3656620  0.57434140  1.03110160
## [5,] -0.357031178 -0.9109974 -0.7144578826  0.9458072 -0.09201943 -0.24641243
## [6,] -0.117014069 -0.2321068  0.8217520978  1.0968101 -1.18647049  1.68628074
##             PC43        PC44        PC45        PC46        PC47       PC48
## [1,] -0.20996256 -0.26979783 -0.02232234  0.19467958 -0.12852110  0.1186386
## [2,]  0.26534725  0.06124778 -0.14509262 -0.20542773  0.74906735  0.6352985
## [3,]  0.25466846  0.34548537  0.58932503 -0.04111156  0.23739990 -0.4620854
## [4,]  0.02436998 -0.72168563  0.51396741  0.85006150 -0.51750603 -0.5766830
## [5,]  1.06357153 -0.47374715 -0.46508358  0.17740669  0.07338641  0.2011894
## [6,]  0.64559540  0.76391268 -0.33056362 -1.40686319  0.23855990  0.5412779
##             PC49        PC50         PC51        PC52         PC53        PC54
## [1,] -0.47408426  0.36632820 -0.221114792  0.08659287  0.202921887 -0.56429022
## [2,]  0.27977737  0.06500505 -0.379707373  0.15280599 -0.158091159  0.10436843
## [3,]  0.31312759 -0.56695086 -0.372736267  1.03517280 -0.003103641 -0.18345652
## [4,] -0.14841188 -0.90026578  0.678598688 -0.42198488 -0.388367929  0.04813490
## [5,]  0.03735946 -0.13589903  0.005542966  0.25711781  0.443969205  0.08099316
## [6,]  0.57913087  0.36442670 -0.210534692  0.08908420 -0.223842488 -0.03031585
##             PC55        PC56       PC57       PC58        PC59        PC60
## [1,]  0.17788012 -0.09741465 -0.1584140 -0.3327349 -0.23543274  0.14137259
## [2,] -0.31194653  0.40355200 -0.8192633 -0.9560810 -0.03005856  0.27905650
## [3,]  0.51685549 -0.01196752  0.8944203  0.1020099  0.19358496  0.01860159
## [4,]  0.24529401 -0.26928696 -0.8567052  0.7879718  0.89692562 -0.15526450
## [5,]  0.07380586 -0.07165184  0.4605282 -0.7086372  0.02123262  0.15349342
## [6,] -0.04679451 -0.61733981  0.3886963  0.1180046  0.15272849 -0.74067126
##            PC61        PC62         PC63        PC64        PC65       PC66
## [1,]  0.5420337 -0.21522110  0.011337783 -0.17855066 -0.03791060  0.2101025
## [2,] -0.2090061  0.28012117  0.407462036 -0.08357081  0.04968771 -0.2213775
## [3,]  0.3609713  0.06828475 -0.003728144 -0.26821311 -0.02681214  0.3520063
## [4,] -0.3382821 -0.63885301  0.708558150 -0.26610003  0.48286393  0.4088681
## [5,]  0.0682367 -0.31584887 -0.551599596  0.10789834 -0.19562580 -0.0128477
## [6,] -0.7194117  0.55382257  0.514415967  0.10134009 -0.30701522 -0.4432713
##             PC67        PC68       PC69        PC70       PC71       PC72
## [1,] -0.51909689  0.28334891 0.09476049  0.19610574  0.1407269 -0.1277971
## [2,]  0.13757464  0.02999191 0.22958571  0.11785784  0.2285787  0.1875041
## [3,] -0.37177927 -0.29498899 0.21763723  0.01009910  0.4531931 -0.4678119
## [4,]  0.05099891  0.32802894 0.52196959 -0.08467014 -0.0582615 -0.2535461
## [5,]  0.04325537 -0.17818549 0.46153578  0.39289479 -0.3359191  0.2128533
## [6,] -0.21770012 -0.21199440 0.19043841  0.02908808 -0.4286533 -0.1956865
##            PC73        PC74        PC75         PC76        PC77          PC78
## [1,]  0.2428455  0.14955193  0.24772571 -0.278433174  0.32258807  0.0009070064
## [2,] -0.1455737  0.01297870 -0.31442006  0.325249630 -0.12592710  0.1187763524
## [3,]  0.1381617  0.02066202 -0.04140399  0.067418905  0.02249998 -0.2817578199
## [4,]  0.1924255 -0.24635625  0.17189179  0.242693556 -0.15010779 -0.5763854872
## [5,] -0.3613821 -0.39889662 -0.15685371 -0.005084836  0.01039949  0.0588177443
## [6,]  0.1170826 -0.10292304 -0.68971469  0.120030215  0.34315601 -0.2274151946
##            PC79        PC80        PC81        PC82         PC83          PC84
## [1,]  0.2072092  0.11045114  0.02290726  0.16893495 -0.106774335  0.0357632989
## [2,]  0.2019984 -0.13363946  0.22413684 -0.03204937  0.183270913  0.0562400806
## [3,] -0.1783240  0.06125581  0.13483458  0.28799251  0.240764564 -0.0007928556
## [4,]  0.1234795 -0.08927775  0.11526411  0.32315921  0.289921213  0.2614749976
## [5,] -0.0677084 -0.12118897  0.36184767 -0.29040697  0.062395112 -0.0646438834
## [6,] -0.3710100 -0.07073274 -0.02863557 -0.23347614 -0.008981002  0.0587340638
##               PC85        PC86        PC87        PC88         PC89
## [1,] -0.1132907008  0.19209582  0.01316851 -0.17345931 -0.006560923
## [2,] -0.1127852532 -0.01852039  0.09825591  0.09264179 -0.041461218
## [3,]  0.0007632933  0.17275546  0.09530554 -0.11461990  0.020477096
## [4,] -0.2483955280 -0.01237575 -0.01099696 -0.11966783 -0.229702465
## [5,] -0.0040571625 -0.44420163  0.02773361 -0.25533816  0.075757828
## [6,] -0.2745976443 -0.21179701  0.13560804 -0.09763305  0.093732040
##              PC90        PC91        PC92        PC93        PC94        PC95
## [1,] -0.003642731 -0.11289932  0.11608906 -0.07704524  0.06614833  0.11492631
## [2,]  0.024193468  0.03816543 -0.11774730  0.03284336  0.22907143 -0.04057139
## [3,]  0.073810100  0.18713160 -0.23797510  0.09322499 -0.10560470 -0.04332543
## [4,] -0.130182703  0.03867649 -0.26777218  0.05498199  0.12958923  0.21716528
## [5,]  0.061212379 -0.14136550  0.04223394  0.03579572 -0.11236475  0.19984461
## [6,] -0.160984487 -0.07197470 -0.02657317  0.15250861  0.15430530  0.08072158
##              PC96          PC97        PC98        PC99       PC100
## [1,] -0.258270737 -0.0362312424  0.03474901  0.08725786 -0.02848885
## [2,] -0.090128855 -0.0537359918 -0.09445251 -0.05327007  0.01088949
## [3,] -0.006743572  0.0120459893 -0.12821049 -0.08188170  0.12048035
## [4,]  0.037614628 -0.0006408959  0.19091908  0.17668717  0.03509228
## [5,] -0.007310014  0.0389636696  0.13666343  0.05765868 -0.10317519
## [6,] -0.088021629  0.1189092184 -0.09280172  0.14622446  0.18399030
##             PC101        PC102       PC103        PC104       PC105       PC106
## [1,] -0.013583776  0.048909193  0.00216378  0.025381950  0.05621436  0.02345814
## [2,] -0.057960400  0.016098193  0.07359925 -0.101541065 -0.06608601 -0.03191740
## [3,]  0.007994631 -0.004963888  0.02902918 -0.076526590 -0.05591075  0.11492641
## [4,]  0.203271380 -0.113055603 -0.16515165 -0.021816263 -0.02176457  0.07787187
## [5,]  0.138790887 -0.085911058  0.04003697 -0.009699707  0.01974313 -0.11142903
## [6,]  0.050447797  0.133785735 -0.10120763  0.155234520 -0.05841889 -0.03275102
##            PC107       PC108       PC109        PC110        PC111       PC112
## [1,] -0.03918637 -0.10011823  0.04367051 -0.078567824  0.061063148 -0.01210633
## [2,]  0.02430551 -0.02982626  0.10875893  0.009612063 -0.075154447 -0.04936937
## [3,] -0.03971283 -0.04659983  0.06379742  0.044902480 -0.053825635  0.04671148
## [4,]  0.13376939  0.01577792  0.20098914 -0.127069247  0.063238208 -0.03424504
## [5,]  0.06745432  0.07870918 -0.03788183 -0.074087792 -0.006017305 -0.02745800
## [6,] -0.08736603  0.17214458 -0.02495514 -0.034747967 -0.002569887  0.09188466
##             PC113       PC114        PC115       PC116        PC117
## [1,]  0.043628685  0.01734119  0.012264508  0.05154895  0.021639315
## [2,]  0.090416102 -0.04063677 -0.051208845  0.02021354 -0.010790308
## [3,]  0.009164668 -0.01981276  0.045332507  0.01182388 -0.085536545
## [4,]  0.018914629 -0.07159628  0.052524427 -0.02034317 -0.012929678
## [5,]  0.009371864  0.02543947 -0.070515577  0.02912714  0.002946295
## [6,] -0.049328896  0.01240519  0.002689908  0.05069294 -0.042221767
##              PC118        PC119         PC120        PC121         PC122
## [1,] -0.0008515265  0.005301082 -0.0004861255  0.005744919 -0.0003524091
## [2,]  0.0411841056  0.007092108 -0.0349009165 -0.070776929  0.0007109495
## [3,]  0.0040223169  0.012616585 -0.0145490756  0.004363147  0.0002528977
## [4,] -0.0304318336 -0.016297073  0.0377117417 -0.009217490  0.0003120880
## [5,]  0.0338064744 -0.030311422 -0.0581924432  0.011809583  0.0006055602
## [6,] -0.0182105313 -0.028032027 -0.0333198643 -0.021873913 -0.0009814977

Let’s revisit the correlation plot. As the principal components are orthogonal, they exhibit no correlation. Consequently, the correlation plot appears perfectly white, except for instances of autocorrelation.

res2 <- cor(pca$x, method = "pearson")
corrplot::corrplot(res2, method= "color", 
                   order = "hclust", tl.pos = 'n')

One challenge of using PCA is deciding how many principal components to retain. Retaining too many may not sufficiently reduce dimensionality, potentially maintaining multicollinearity. Conversely, keeping too few can lead to significant information loss, reducing the accuracy of your model. A widely used method to determine the optimal number of principal components is the scree plot. This plot displays the eigenvalues of the principal components in descending order and identifies where an abrupt change in the slope occurs, signaling a drop in the importance of subsequent components. The ideal pattern is a steep curve, followed by a bend, and then a straight line. Use the components in the steep curve before the first point that starts the line trend.

factoextra::fviz_eig(pca, barfill = "gray73", barcolor = "gray73",
                     linecolor = "steelblue2", xlab = "Principle Component", ncp = 30)

As shown in the scree plot above, at first glance, it appears that the curve starts to level off from PC10, suggesting that including the first 9 components may be optimal. According to the summary results of the PCA object, including these first 9 components explains 70.85% of the variance.

Now let’s rerun the model with 9 principal components and compute the VIF scores for each component. As shown below, all the VIF scores are 1, confirming the absence of multicollinearity as previously indicated by the correlation matrix. Instead of including numerous correlated variables, which could lead to inefficient estimates, the revised model now offers a more efficient set of estimates. Moreover, it still explains a significant portion of the variance in the outcome variable, with an \(R^2\) of 0.6261, only slightly lower than the 0.6979 of the original model.

crimeData.pca <- as.data.frame(cbind(y = crimeData.y, pca$x[, 1:9]))
fit_pca <- lm(y ~ ., data = crimeData.pca)
summary(fit_pca)$r.squared

## [1] 0.6260937

round(vif(fit_pca), 2)

## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 
##   1   1   1   1   1   1   1   1   1

For a more stringent method to determine how many principal components should be included in the linear regression model, please see here for an example.

Another challenge of using PCA is to interpret the meaning of the principal components. Unlike the original features, which may have a clear definition and unit, the principal components are abstract and hard to explain. One way to interpret the principal components is to look at the coefficients of the linear combination that forms each component. These coefficients indicate how much each original feature contributes to the component and what is the direction of the relationship. For example, if a component has a high positive coefficient for a feature, it means that the component increases with the feature.

The function fviz_contrib function from factoextra package can be used to draw a bar plot of variable contributions. If your data contains many variables, you can decide to show only the top contributing variables. The R code below shows the top 50 variables contributing to the first principal component we have in this example:

factoextra::fviz_contrib(pca, choice = "var", axes = 1, top = 50, 
                         fill = "gray73", color = "gray73")

The red dashed line on the graph above indicates the expected average contribution. If the contribution of the variables were uniform, the expected value would be 1/length(variables) = 1/122 = 0.82%. For a given component, a variable with a contribution larger than this cutoff could be considered as important in contributing to the component.

For instance, in the linear regression model where we include principal components instead of the original variables, there is a negative relationship between our outcome variable, violent crimes per capita, and PC1. According to the contribution plot, variables such as x.V25, x.V18, x.V50, and others above the red dashed line indicate that the relationship between our outcome variable and these variables is likely negative.

Discussion V

TBD

Discussion VI

TBD